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The long way of trying to reproduce on the Earth the processes which make the stars brighten, 
is the more than 50 years long history of nuclear fusion. Despite the last half century has 
seen an extraordinary improvement of the scientific knowledge and the associated human 
technological capabilities, the goal of employing the nuclear fusion reactions as profitable 
and peaceful energy source, has still to be achieved. This large amount of efforts and time 
that have been spent may suggest that achieving nuclear fusion on the Earth is not an easy 
task. Indeed, it is not. 

Fusion relics on the nuclear interaction between light elements which are combined to 
form heavier ones, releasing an amount of energy. On the other hand, such interactions are 
possible if the original components have an high free energy, thus implying that the reacting 
medium is locally far from the thermodynamical equilibrium. Practically, that means that 
the mixture of the reacting components, typically deuterium and tritium, have enough energy 
to be completely ionized, thus forming a plasma: external constraints have then to balance 
its internal pressure. At the very beginning of this work, it is important to highlight that 
the actual crucial issue of the nuclear fusion can be extremely summarized by a single word: 
confinement. Prom the problem of the energy confinement in fact, stems a wide galaxy of 
both scientific and technological topics which are subjects of the living history of the nuclear 
fusion research. Even if this thesis will only deal with a very specific issue in this long list, the 
main concern that has inspired this work was to clearly locate the research activity within 
this global framework. Understanding, predicting and controlling the conGncment mecha- 
nisms are the main keys to thrive over the scientific and technological challenges of fusion. 

A bit more quantitatively, a characteristic time scale can be defined as the ratio between 
the energetic content of the plasma and the power losses: this defines the so called energy 
confinement time te- In the case of the tokamak configuration, where a particular mag- 
netic topology of nested flux surfaces provides a pressure that counterbalances the internal 
plasma one, the losses are mainly due to conductive effects. The energy confinement time can 
then be associated to a diffusive time a/xth) where a is the typical macroscopic scale of the 
plasma column and Xth is the effective thermal diffusivity. Leaving aside any technological 
implication, which de facto represents a wide and very active domain of research within the 
fusion community, this study will treat the issue of confinement only in terms of a problem 
of transport of energy and matter. 

The transport in tokamak plasmas is usually referred as anomalous. The anomaly makes 
reference to the experimental evidence, verified on many devices since the 1980s, that the 
measured thermal diffusivities are larger by one or two order of magnitudes, respectively 
for the ions and the electrons, than the expectations based only on the coUisional processes 
(neoclassical theory). Effectively, the anomalous tokamak transport is largely due to the 
micro-turbulence that affects the plasma dynamics. The main drives of these turbulent 
micro-instabilities were already known in the 1960s- 70s: these are namely the temperature, 
density and velocity gradients transverse with respect to the magnetic flux surfaces. The 
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resulting instabilities originate fluctuations in the plasma pressure as well in the electromag- 
netic field associated to the spatial charge distribution. The final macroscopic result is the 
presence of turbulent fluxes of particle, energy and momentum. 

The quest for comprehensive and reliable models of the tokamak anomalous transport 
lasts from several decades. From one side, the combined efforts of sharing the increasingly 
large amount of experimental results collected from many fusion facilities across the world, 
allowed the formulation of semi-empirical scaling laws. The underlying idea of this approach 
is to try to identify a statistically coherent correlation, which is therefore transferred into 
a predictive capability, between a number of physically relevant dimensionless parameters 
(related for ex. to the size of the device, the plasma pressure, etc.) and the critical fusion 
quantities, typically the energy confinement time. On the opposite side, the ambition of the 
flrst principle modelling is to elaborate a theoretically based understanding of the plasma 
dynamics. Because of the complexity of the underlying processes, which often makes that the 
problem cannot be analytically solved, the first principle quantitative information is more 
and more deferred to the numerical simulations, strongly pushed by the recent advances of 
the computational capabilities. These two different approaches, the experimentally driven 
formulation of scaling laws and the coupled theoretical-numerical effort for first principle 
models, refer respectively to an inductive rather than a deductive knowledge process, or in 
other words, to a top-down rather than a bottom-up strategy. Not pretending to cpiste- 
mology, the scientific understanding does not rely on a purely inductive neither deductive 
reasoning, but on a balanced hypothetical-deductive method. The work here proposed aims 
developing a reduced transport model for the core of tokamak plasmas which stems from the 
combined contributions of theoretical, numerical and experimental insights. 

A rigorous first principle modelling of the tokamak turbulence should deal with the self- 
consistent dynamics between the electromagnetic fields and the plasma response with respect 
to the field perturbations. The resulting problem is of great complexity, since it involves a 
nonlinear dynamics in both the spatial and the velocity phase spaces. The actual state of 
the art is represented by the coupling of the Maxwell equations with a reduced 5-dimensional 
description of the plasma response: this is the so called gyro-kinetic approximation, which 
in fact results from the average over the fast particle gyro-motion. Moreover, a number of 
additional simplifications arc usually adopted, provided the fact that the plasma dynamics 
is intrinsically extended over disparate temporal and spatial scales. Still, the quantitative 
information that can be gained from the numerical solution of such an electromagnetic non- 
linear system requires huge computational resources and it is hardly applicable for inferring 
predictions over the experimental macroscopic time scales of a typical tokamak discharge. 

The approach of this thesis work adopts a different strategy. The turbulent transport 
model here proposed follows from a gyro-kinetic quasi-linear approximation. The relevance 
of such a solution goes well beyond the motivations of simplifying the nonlinear nature of the 
problem and therefore lowering the computational cost: the reduction of complexity under 
a minimal set of hypotheses is a physically challenging issue. This logic of modeling allows 
identifying whether the tokamak turbulence is reminiscent or not of the linear dynamics, 
resulting in a great progress of the physical understanding. 

A linear response that keeps the fundamental kinetic wavc-particlc resonance mechanism 
in the velocity space is chosen. On the other hand, the nonlinear self-consistency of the 
system is broken and split into two basic parts: (1) the linear response of the plasma and 
(2) the saturated fluctuating potential. The hypotheses underlying the linear response, and 
therefore its limits, have to be accurately verified against fully nonlinear simulations. The 
saturation regime can instead be studied by mean of both turbulence measurements in the 
plasma core, and nonlinear simulations. At any level of the whole process of the model for- 
mulation, simplified analytical models can greatly enhance our understanding of both the 
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experimental, hence real, and the numerical, hence artificial, observations. 
The final goal of such a reduced quasi-linear gyro-kinetic model is the solution of the time- 
dependent transport problem. Finally, the detail of the energy and particle turbulent trans- 
port at the microscopic scales predicted by the model, is recast into an upper-level solver 
which integrates all the information concerning the equilibrium, the sources and the transport 
of the plasma: the time evolution of the thermodynamically relevant quantities is therefore 
predicted. 

In conclusion, this thesis work tackles the issue of the energy and particle confinement in 
the core of tokamak plasmas through the formulation of a reduced physical model: the related 
understanding, validity and limits stem from the concurrence of the theoretical, experimental 
and numerical analysis. 
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Le long chemin d'essayer de reproduire sur Terre les processus qui font briller les etoiles, 

est I'histoirc, longuc de plus que 50 ans, de la fusion nuclcaire. Bien que la dernicre moitie 
du siecle ait connu une extraordinaire amelioration de la connaissance scientifique et des ca- 
pacites technologiques associees, I'objectif d'utiliser les reactions de fusion nuclaire en tant 
que soiucc d'cncrgic rentable et pacifique, n'a pas encore etc atteint. Cettc grandc quantitc 
d'efforts et de temps qui a ete depensee suggere que la realisation de la fusion nucleaire sur 
Terre n'est pas une tache facile. 

La fusion repose sur I'interaction a rcchcUc nuclcaire cntre des ions Icgers qui se com- 
binent pour former des ions plus lourds, en liberant de I'energie. D'autre part, cette fusion 
est possible si les composants d'origine ont une energie libre elevee, ce qui implique que ce 
moyen rcactif est au niveau local loin de I'cquilibrc thermodynamique. En pratique, cela 
signifie que les composants qui interviennent dans les reactions, generalement le deuterium 
et le tritium, ont sufiisamment d'energie pour etre completement ionisees, formant ainsi un 
plasma: des contraintes cxterncs doivent done pourvoir a cquilibrer la prcssion interne. II est 
important de souligner que la question cruciale de la fusion nucleaire pent etre synthetisee 
par le mot conSnement. A partir de la question du confinement de I'energie en fait, resulte 
une grande variete de problcmes a la fois scicntifiques et technologiques qui constitue les 
sujets de I'histoire de la recherche sur la fusion nucleaire. Meme si cette these ne portera que 
sur une question tres precise dans cette longue liste, la principale preoccupation qui a inspire 
ce travail a ete de bien localiser I'activitc de recherche dans ce cadre global. Comprcndre, 
prevoir et controler les mecanismes du confinement sont les cles principales pour reussir dans 
les defis scientifiques et technologiques de la fusion. 

De facon un pen plus quantitative, une cchelle de temps caractcristique du confinement 
peut etre definie comme le rapport entre le contenu energetique du plasma et la puissance 
perdue: cela est la definition du temps de confinement de I'energie te- Dans le cas de la 
configuration tokamak, oil une topologie magnetique des surfaces de flux imbriquees four- 
nit une pression qui contrebalance celle interieure du plasma, les pertes dans le coeur sont 
principalement dues a des effets de conduction. Le temps de confinement de I'energie peut 
alors etre associe a un temps de diffusion a/xth, ori a est la grandeur macroscopique de la 
colonne de plasma et Xth est la diffusivite thermique efficace. Laissant de cote toute impli- 
cation technologique, qui represente de facto un domaine large et actif de la recherche dans 
la fusion, cette etude traite la question du confinement seulement en termes d'un probleme 
de transport de I'energie et de la matiere. 

Le transport dans les plasmas de tokamak est generalement qualifie comme anormal. 
Cette anomalie fait reference a I'evidence, verifiee sur des nombreuses installations depuis 
les annees 1980, que les diffusivites thermiques sont plus grandes de un ou deux ordres de 
grandeur, respectivement pour les ions et les electrons, des predictions basees uniquement 
sur les processus de collision (theorie neoclassique). En effet, le transport anormal dans les 
tokamaks est largement du a la micro-turbulence qui affecte la dynamique du plasma. Les 
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mccanismcs principaux a la base dc ccs micro-instabilits ctaicnt dcja connus dans Ics annccs 
1960-70: ce sont notamment les gradients de temperature, densite et vitesse perpendiculaires 
aux surfaces de flux magnetique. Ces instabilites sont a I'origine de fluctuations de la pres- 
sion du plasma ct du champ clcctromagnctiquc associc a la distribution spatialc de charge. 
Le resultat macroscopique est la presence des flux turbulents de particules, d'energie et de 
moment angulaire. 

La quote dc modclcs cxhaustifs ct fiablcs du transport anormal dans Ics tokamaks dure 
depuis plusieurs decennies. D'un cote, les efforts combines de partage d'une grande quantite 
de resultats experimentaux obtenus aupres des nombreuses installations de fusion a travers 
Ic mondc, a pcrmis I'claboration des lois d'cchcllc scmi-cmpiriqucs. Ccttc approchc pcrmct 
de trouver une correlation statistique coherente, qui est transferee en capacite de prediction, 
entre un certain nombre des parametres physiques adimensionnels (lies par ex. a la taille de 
I'apparcil, la prcssion du plasma, etc.) ct des quantites critiques pour la fusion, gcneralement 
le temps de confinement d'energie. De I'autre cote, I'ambition de la modelisation de la turbu- 
lence consiste a acquerir une comprehension de la dynamique du plasma basee sur la theorie. 
En raison dc la complcxitc des processus non-linoairos sous-jaccnts, qui fait souvcnt que Ic 
probleme ne pent pas etre resolu analytiquement, seule la simulation numerique, fortement 
poussee par les recents progres des capacites de calcul, permet de donner des reponses quan- 
titativcs. Ccs deux approchcs, la formulation des lois d'cchcllc et Ics efforts thcoriqucs ct 
numeriques, sont respectivement des processus de connaissance inductive et deductive, soit 
des strategies top-down plutot que bottom-up. Sans pretention d'epistemologie, la connais- 
sance scicntifiqiie ne sc fondc pas sur un raisonnement purcmcnt inductif ni dcdiictif, mais sur 
une methode hypothetique-deductive. Le travail propose ici vise a I'elaboration d'un modele 
de transport reduit pour le coeur des plasmas de tokamak, qui integre des contributions a 
caractcrc thcoriquc, numerique ct experimental. 

Une rigoureuse modelisation aux premiers principes de la turbulence dans les tokamaks de- 
vrait traiter la dynamique auto-consistante entre les champs electromagnetiques et la reponse 
du plasma a I'egard des perturbations du champ. Le probleme qui en resulte est d'une grande 
complexite, car il implique une dynamique non-lineaire a la fois dans les espaces de phase 
en espace et en vitesse. L'etat de Fart actuel est represente par le couplage des equations 
de Maxwell avec une description rcduitc a 5 dimensions de la reponse du plasma: c'est 
I'approximation dite gyro-cinetique, qui resulte de la moyenne sur le mouvement de gira- 
tion rapide des particules le long du champ magnetique. En outre, un certain nombre de 
simplifications supplcmentaires est gcneralement adoptc, utilisant le fait que la dynamique 
du plasma est intrinsequement etendue sur des echelles temporelles et spatiales disparates. 
Toutefois, I'information quantitative qui peut etre retiree par la solution numerique d'un 
tel systeme clcctromagnctiquc non-lineaire exige d'enormcs ressources de calcul, ct elle est 
actuellement difficilement applicable pour obtenir des predictions sur les echelles de temps 
macroscopiques d'une decharge typique de tokamak. 

L'approche dc ce travail dc these adoptc une strategic differente. Le modele de transport 
turbulent ici propose part d'une approximation gyro-cinetique quasi-lineaire. L'importance 
d'une telle solution va au-dela des raisons de simplification de la nature non-lineaire du 
probleme afin de reduire le cout de calcul. En cffct, ccttc logique dc modelisation permet 
d'identifier dans quelle mesure la turbulence des tokamaks montre des reminiscences de la 
dynamique lineaire, entrainant une progression sensible de la comprehension physique. 

Une reponse lineaire qui conserve Ic mecanisme fondamental de resonance cinetique 
onde-particule dans I'espace des vitesses est choisie. D'autre part, I'auto-consistance non- 
lineaire du systme, ici defectueuse, est scindee en deux termes : (1) la reponse lineaire du 
plasma et (2) le potentiel electrostatique sature fluctuant. Les hypotheses qui sous-tendent la 
reponse lineaire, et done ses limites, doivent etre attentivement confrontees a des simulations 
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entierement non-lincaires. Le regime dc saturation peut ctrc ctudie a la fois par moycnne 
des mesures de turbulence dans le coeur du plasma, et par des simulations non-lineaires. 
A tout niveau de I'ensemble du processus de la formulation du modele, des modeles analy- 
tiqucs simplifies pcuvent grandcment amcliorcr notrc comprehension dc la complexite a la 
fois experiment ale, done reelle, et numerique, done artificielle, des observations. 

En conclusion, ce travail de these aborde la question du confinement de I'energie et des par- 
ticulcs dans le coeur des plasmas dc tokamak a travers la formulation d'un modele physique 
reduit: la validite et les limites de ce modele decoulent de la concurrence entre I'analyse 
theorique, experimental et numerique. 

Contenu de la these 

Le contenu de ce travail de these est le suivant. 

Dans le Chapitre 2, la strategic de confinement adoptee par les tokamaks est introduite. 
Les principales instabilites du plasma, responsables du transport turbulent de I'energie et 
la maticrc dans un tcl systcmc, sont traitccs. Les deux representations fondamentalcs du 
plasma, celle fluide et celle cinetique, sont ici brievement presentees. Une attention parti- 
culiere est consacree aux raisons pour lesquelles une approche gyro-cinetique a ete preferee 
dans la modclisation quasi-lincairc. Un exemple pertinent pour le tokamak est prcscntc afin 
de souligner I'importance de retenir correctement la resonance cinetique onde-particule: il 
s'agit de la dependance du seuil d'instabilite lineaire en fonction du rapport des temperatures 
Tj/Te des modes clcctroniques et ioniques. 

Le Chapitre 3 traite la question de la reponse quasi-lineaire. Premierement, la derivation 
du modele, appele QuaLiKiz, et ses hypotheses sous-jacentes permettant d'obtenir des flux 
turbulents d'cncrgic et dc particulcs, sont presentees. Deuxicmemcnt, la validite dc la reponse 
quasi-lineaire est confrontee a des simulations gyro-cinetiques non-lineaires afin de: (a) iden- 
tifier les temps caracteristiques dominant la dynamique turbulente, (b) comparer les relations 
de phase entre les champs fluctuants, (c) examiner I'intensite globale du transport, normalise 
a I'intensite du potentiel sature. 

Le Chapitre 4 porte sur le modele de la saturation non-lineaire. Dans la premiere partie, 
les simulations gyro-cinetiques non-lincaires sont validees quantitativement par rapport aux 
mesures de turbulence dans le tokamak Tore Supra. Les spectres des fluctuations de densite 
tant dans I'espace, k, que en frequence, oj, sont etudies aussi en termes de modeles analy- 
tiqucs simplifies. Sur la base de cctte validation non-lineaire avec les mesures, le modele de 
saturation qui est introduit dans QuaLiKiz est presente et discute. 

Le Chapitre 5 est consacre a qualifier les resultats de QuaLiKiz. Les flux quasi-lineaires 
d'cncrgic et des particulcs sont compares aux predictions des simulations non-lincaires pour 
une large gamme des parametres. Finalement, le couplage de QuaLiKiz dans le solveur de 
transport integre CRONOS est presente. Cette procedure permet de resoudre la dependance 
temporelle dans le probleme de transport, et done une application directe du modele a 
I'experience. Des resultats preliminaires concernant I'analyse experimentale sont enfin ex- 
amines. 
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1.1 Basics around nuclear fusion 

The nuclear fusion reactions rely on the interaction between light nuclei on the scale of 
characteristic range of the nuclear strong force: the overcome of the electrostatic Coulomb 
barrier gives rise to the fusion of the nuclei to form heavier elements and releasing the 
difference of the binding energy of the original components. Because of the presence of the 
Coulomb repulsion, the nuclei involved in the fusion reaction must have a kinetic energy of 
the order of energy barrier: the elements will therefore form a completely ionized plasma. 

The nuclear fusion processes are extremely common in the universe as they represent 
the source of energy of the stars. Inside the sun, for example, a series of nuclear reactions 
converts mainly hydrogen nuclei into helium ones through the so called proton cycle, whose 
net result is: 

4p+ =^ ^iJe + 2e+ + 2ty + 27 MeV (1.1) 

These are processes that have a particularly high energy efficiency, but the drawback of an 
extremely low reaction rate. This is not a problem in the stars, because the extreme pressure 
of the core provides a gravitational confinement of the particles for a sufficiently long time, 
so that a balance is established between the radiation pressure due to the thermonuclear 
reactions and the gravitational compression. 

In order to reproduce on the earth controlled nuclear fusion reactions, it is necessary to 
adopt appropriate solutions both in terms of the nuclear reactivity and the plasma confine- 
ment. Although from the theoretical point of view several fusion reactions can be explored, 
as shown in Fig. |1.1[ the way that nowadays appears as the most accessible and promising 
one, relies on the fusion between between deuterium D and tritium T nuclei: 

^D+^T => n + '^i/e + 17.6 MeV (1.2) 

As shown in Fig. |1.1[ the D-T reaction maximizes the reactivity at the lowest ion temperature, 
with a cross-section which is on average higher by two orders of magnitude compared to the 
other reaction channels. It is interesting to note that for the most part of the elementary 
fusion reaction processes, the overcome of the Coulomb barrier between the nuclei occurs 
through a quantum tunneling effect. 

The D-T reactions are also interesting when considering the coupling with the following 
one: 



n+ ^Li 



+ ^He + 4.8 MeV 



(1.3) 
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Figure 1.1: Reactivity (reaction cross-section) for different nuclear fusion processes. 



The nuclear reactions |1.1| and |1.2| suggest in fact that the energetic neutrons produced by 
D-T reaction |1.1[ could interact, within an appropriate structure that is usually referred to 
as blanket, with lithium atoms, thus regenerating the tritium. This would end in a twofold 
advantage: the first one relates to the possibility of adopting a closed cycle for the tritium in 
a future fusion reactor, while the second one refers to the creation of a thermal energy source 
associated with the neutron-heated lithium, hence usable for electricity production purposes. 
The internal total energy of a thermonuclear plasma can be defined as: 



W 



The famous energy confinement time can then be formally expressed as: 



p 



w 

- dW/dt 



(1.4) 



(1.5) 



where Pin indicates the external input power. The break-even condition simply foresees that 
the nuclear fusion power overcomes the external input power: this is usually expressed by 
the so-called Q-factor, i.e. Q — ^-p^ > 1. On the other hand, the ignition corresponds to a 
self-sustaining reaction, i.e. Q — > cxi, which physically means that the nuclear fusion power 
associated to the a-particles balances the power losses. Historically, this criterion has been 
translated into an operational condition, the famous Lawson's triple product: 



n,T,Te > 3 • 10^^ m"3 keVs 



(1.6) 



The ion temperatures which are required are typically of the order of w 10 — 30 keV. 
Seeking the maximization of the product of the plasma density and the energy confinement 
time leads to two different strategies, that are at the present day the two ways the scientific 
community follows for achieving the nuclear fusion on the Earth. 

The way of the inertial confinement plans to maximize the ion density, at the detriment 
of the energy confinement time. Practically, the inertial fusion most probably consists of a 
pulsed regime of micro-esplosion of small D-T targets; the ignition is triggered, according to 
several different schemes, by compressional waves induced with powerful lasers or accelerated 
particle beams. 

The second approach is the magnetic fusion. In this case, the ion density is much smaller 
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with respect to a gas at the atmospheric pressure, but the energy confinement time aims to 
reach macroscopic scales of the order of the s. An appropriate magnetic configuration, which 
is not necessarily the tokamak one, is responsible to spatially confine the hot temperature 
plasma. 

1.2 Outline of this work 

The outline of this thesis work is the following. In Chapter [2] the framework of the tokamak 
strategy to deal with confinement, hence the main plasma instabilities which are responsible 
for turbulent transport of energy and matter in such a system are introduced. The two 
principal plasma representations, the fluid and the kinetic ones, are here briefly introduced. 
Particular attention is devoted to the reasons why a gyro-kinetic approach has been preferred 
in the quasi-linear modelling. A tokamak relevant case is presented in order to highlight 
the relevance of a correct accounting of the kinetic wave-particle resonance: the example 
deals with the T^/Tg dependence of the linear instability threshold for the ion and electron 
modes. The Chapter |3] discusses the issue of the quasi-linear response. Firstly, the derivation 
of the model, called QuaLiKiz, and its underlying hypotheses to get the energy and the 
particle turbulent flux are presented. Secondly, the validity of the quasi-linear response is 
verified against the nonlinear gyrokinetic simulations in order to: (a) identify the dominant 
characteristic times of the turbulent dynamics, (b) compare the phase relations between 
the fluctuating fields between the linear and the nonlinear phase, (c) examine the overall 
transport intensity, normalized to the intensity of the saturated potential, between the linear 
and the nonlinear regime. The Chapter |4] deals with the model of the nonlinear saturation. 
In the first part, the nonlinear gyrokinetic simulations are quantitatively validated against 
turbulence measurements on the Tore Supra tokamak. Both the spatial k and the frequency 
Lu spectra of the density fluctuations are investigated also in terms of simplified analytical 
models. Consequently, the saturation model that is assumed in QuaLiKiz is presented and 
discussed. The Chapter [5] is devoted to qualify the global outcomes of QuaLiKiz. Both the 
quasi-linear energy and the particle flux are compared to the expectations from the nonlinear 
simulations, across a wide scan of tokamak relevant parameters. Therefore, the coupling of 
QuaLiKiz within the integrated transport solver CRONOS is presented: this procedure allows 
to solve the time-dependent transport problem, hence the direct application of the model to 
the experiment. The first preliminary results regarding the experimental analysis are finally 
discussed. 
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Chapter 2 

Description of the plasma 
instabilities 

2.1 General framework 

The purpose of this first section is to give a general overview of the main instabihties which 
can drive anomalous transport in magnetically confined plasmas. This framework is intended 
to introduce the basic mechanisms which are modeled in this thesis work and will be exten- 
sively discussed in the following sections. Even if specific attention will be devoted to the 
toroidal magnetic configurations, this description should not suffer of a loss of generality. In 
fact, the two main instabilities detailed here below, namely the interchange and drift-waves, 
can commonly characterize the turbulence of electrically charged media in the presence of a 
eventual background magnetic field. 

In the first paragraph, a brief description of the charged particles trajectories in toroidal 
magnetic geometries is given; their characteristic drift velocities and the invariants of motion 
are introduced. 

In the second paragraph the main physical features of the interchange and drift-wave instabil- 
ities are presented. The first mechanism is formally analogous the hydrodynamic Rayleigh- 
Bernard instability, and can take place even in two dimensions thanks to the non-homogeneity 
of the magnetic field. The second one arises from both the single particle drifts and the col- 
lective behavior of the plasma. This is essentially a three dimensional effect, relying on the 
plasma response to the perturbations in the direction parallel to the magnetic field. 



2.1.1 Geometry and particle motion 

The tokamak magnetic geometry is axis-symmetric and consists of a series of closed nested 
surfaces. The toroidal component of the magnetic field is produced by external coils, while 
the poloidal field is originated by the current induced in the plasma. The resulting magnetic 
field is usually expressed by the relation: 

B = /(V')V(/7 -f VV- X (2.1) 



In Eq. (2.1 1 ■0 is the magnetic poloidal flux normalized to 27r, while ip stands for the toroidal 
angle and /(V') is a flux function. The field lines are winded around the magnetic flux surfaces. 
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An important parameter for tokamak plasmas defines the winding rate of these field lines, 
constant on a given flux surface, in the following way: 

where 9 is the poloidal angle. g(V') is the so called safety factor. Introducing the radial co- 
ordinate r that labels a given magnetic flux surface through a simple dependence of the flux 
function / = /(r), all the quantities of interest can the be decomposed on the orthonormal 
basis (er,eg,e<^). 

The most recent tokamak configurations adopt a plasma cross-section that presents fi- 
nite elongation and triangularity, with the typical D shaping, leading to improved stability. 
Nevertheless, a more simple circular geometry is still of great interest. Within this approxi- 



mation, exemplified in Figure 2.1 the flux surfaces are assumed circular and concentric, and 
the magnetic field can be simply written as B = Bgeg + B^e^p, giving: 



where the is the reference field value on the magnetic axis and i?o is the major radius 
of the plasma. In this context, it is worth noting that the magnetic field B can be as well 




Figure 2.1: Scheme of the simplified magnetic geometries, characterized by circular concentric 
flux surfaces 



expressed through the so called Clebsch representation [65]: the latter is characterized by 
the definition of the new variable C, = — q{ip)9. This advantage of this procedure follows 
from the fact that the magnetic field can now be written as: 

B = V?7 X VtA (2.4) 

which satisfies B • VC — B • Wip — 0. The new coordinates system (-0,6', C) appears then to 
be a field-aligned representation, where the variable 9 refers simultaneously to (a) an angle 
in the poloidal plane (at fixed ip), or (b) a parametrization of distance along the field line (at 
fixed C.). The field-aligned coordinates systems appears particularly useful when applied to 
numerical simulations of tokamak microturbulence: examples will be provided in the follow- 
ing of this thesis work. 
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Particle motion 

Under the effect of the just described magnetic field, a charged particle with mass m 
and electric charge e undergoes a cyclotron gyromotion, which is characterized by typical 
spatial and temporal scales: these are respectively the Larmor cyclotron gyroradius pc = 
mv±/eB and the gyrofrequency ujc = eB/m, where v±_ is the amplitude of the particle 
velocity transverse to B. The motion of this charged particle in a arbitrary time-dependent 
electromagnetic field cannot be analytically solved. On the other hand, these trajectories are 
found to be integrable under the hypothesis of scale separation. The key concept relies on 
the decoupling of the fast cyclotron dynamics in the plane transverse to B, from the slow 
motion of the particle guiding center. In magnetic fusion applications, this limit appears 
to be a very good approximation, since the external field B can be assumed quasi-static in 
terms of spatial and time variationf[^ In particular, the following ordering is established by 
the adiabatic theory [75] : 

• The typical spatial variations of the magnetic field are negligible compared to the 
particle gyroradius 

VB , , 

(2.5) 

• The typical time variation of the magnetic field can be neglect in comparison to the 
particle gyrofrequency 

^«v.VB«^i?«c.,S (2.6) 
at rto 

where Vty^ is the thermal speed of the particle. 

Within this framework, i.e. with the hypothesis of a quasi-static electromagnetic field, 
three fundamental adiabatic invariants are conserved along the particles trajectories; these 
are: 

• The total energy of the particle E 

E = ^mwg|| + fiB + e(j) (2.7) 

where the presence of (j) refers to the eventual presence of stationary electric potential 
and the index G recall the guiding center framework. 

• The magnetic moment fi 

In other words fj, expresses the magnetic flux across the particle gyromotion; a direct 
consequence of the invariance of this quantity is the variation of the Larmor gyroradius 
along the field line. 



The kinetic toroidal moment M 



M = e^ + mRvG4, (2-9) 



^An important exception that will not be here treated is related to the plasma heating through RF power 
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It is possible to show that this invariant is a direct consequence of the axis-symmetry of 
the system. Within the Hamiltonian mechanics formulation M is in fact the conjugate 
momentum of the toroidal angle if, M = dL/dip, where the Lagrangian of the particle 
L — l/2mv^ — e(f> — eA • v (A is the potential vector) does not depend on (p. This 
is an idealization of the real tokamak plasmas, where inhomogeneities in the toroidal 
magnetic field are present due to the finite numbers of coils. 

At this point it is particularly useful to use the just described invariants in order to introduce 
the normalized velocity-space coordinates, defined as 

£=^ (2.10) 

(2.11) 

^ = sgn(«,j|) (2.12) 

such that (f , A) are again unperturbed constants of motion, while c, is the sign of the parallel 
velocity. This notation helps in highlighting a general feature of the dynamics of charged 
particles in non-uniform magnetic fields, i.e. the magnetic mirror phenomenon. Due to the 
radial dependence of B in tokamak plasmas, some particles have not enough energy in the 
parallel direction to undergo a complete turn in the poloidal direction. Writing the parallel 
velocity for a given particle with the new notation, we have 



"II 



J—y/£^l-\h{r, 9)^ (2.13) 

V m-s 



where b{r,9) = B{r,9)/B{r,0). It appears that the particles satisfying the condition A = 
b{r,9o)~^ have a bouncing motion along the field line between the poloidal angles 9 = —9o 
and 9 ~ +9q. Two main classes can then be distinguished, passing and trapped particles. 

The fundamental equation governing the guiding center dynamics can be derived from 
the non-relativistic Lorentz equation of motion: mdv/dt = e(E -f v x B). The latter can be 
adopted to the evolution of the guiding center with an extra term which embeds the effect 
of the fluctuations of B at the Larmor scale, giving: 

dvr' 

TO-^ = e(E + VG X B) - /iVB (2.14) 



The projection of the Eq. (2.14) on the transverse and parallel directions allows the under- 
standing of some essential features of the particles dynamics in tokamak plasmas. Starting 
with the expansion on the perpendicular plane, four basic drifts mechanisms for the guiding 
center motion are identified: 



The E X B drift 



E X B , , 

VExB = ^2 (2-1^) 



• The VB drift 
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The curvature drift 



B 



R 



where 



N 



V±B 
B 



(2.17) 



Introducing the dimensionless ratio between the kinetic and magnetic plasma pressure 
P = p/ usually of the order of 0.01 in tokamaks, the second term of the 

quantity N/i? is negligible, so that N/R w V^_B/B. 



The polarization drift (higher order) 



toB 
eB2 



(vexB + VvB + Vc 



(2.18) 



The appearance of these drift motions define some fundamental properties of the tokamak 
plasma instabilities described in the next paragraph. 



2.1.2 Drift waves and interchange instabilities 

The purpose of this paragraph is only to show the main physical mechanisms at the origin of 
the turbulence object of the modeling presented in this work. In tokamak plasmas there are 
in fact generally many free energy sources, originating a very wide spectrum of micro/macro- 
instabilitics. whose exhaustive discussion is far beyond the scope of the present work. A 
general distinction can be made on the pressure-driven modes or the current-driven modes. 
The latter are typically driven by the current flowing in the parallel direction and they are 
usually described within the framework of fluid MHD [Magneto Hydro Dynamics) models. 
Even if this class of instabilities has important consequences in tokamak plasmas, setting for 
example intrinsic limits on the total plasma current and pressure, these modes will not be 
treated here. Indeed, the quasi-linear turbulence modeling presented in this work presupposes 
a flxed magnetic equilibrium, where the pressure driven micro-instabilities are completely 
decoupled from the evolution of the macroscopic magnetic fleld. The self-consistent interplay 
between the MHD and the pressure driven instabilities has just started to be explored, and 
represents one of the next research challenges in the progress of understanding magnetic 
fusion. 

The drift-wave instabilities arc a key mechanism in tokamak plasma turbulence. In the 
case of electrostatic turbulence, the transport is set by the fluctuations in the E x B drift 
velocity 

*v. = ^ (2.19) 

The thermal velocities of electrons differs from the ions one by a factor y^rni/m^; in first ap- 
proximation then, the electrons can be supposed to instantaneously respond to the potential 
fluctuations characterized by a smaller frequency with respect to the their parallel dynamics, 
i.e. oj <C k\\Vth,e- In the parallel force balance equation for electrons 

dVp B femeiw 

mene— ^ • 15 = e"<.V||0 - V||Pe + 2.20 
dt B e 

the dominant terms arc the electric field and the pressure gradient. This relation reduces 
then to V|ipe = e«eV|i0, when neglecting the temperature fluctuations and in the collisionless 
limit. A linear Boltzmann response directly follows as: 
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Eq. (2.21) is usually referred as the hypothesis of adiabatic electron response. Following a 

(2.22) 



fluid approach fl4j, the fluid dynamics equation is: 

mn (St + u • V) u = nq (-V0 + uxB)-Vp-V-n 



where u and v are respectively the fluid an^ kinetic velocity while the tensor 11 contains the 
non-diagonal terms of the pressure tensor P = pi + 11 ^ ^ J d?v (v — u) ® (v — u) /. At the 
first order of the expansion in the parameter e = — ^ u±k±_ ^ ^j^^ perpendicular projection 



of the Eq. \2.22\ reduces to nq (— V0 + u_l x B) — V_lp — 0. Hence, the perpendicular fluid 

B X V0 



velocity at the first order reads: 



Ufi + u = 



B Vp 

en 



(2.23) 



where u* is called the diamagnetic drift velocity. 
Under the hypothesis of adiabatic electron response described by Eq. (2.211, the electric 

(a) Drift wave (b) Drift instability 
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Figure 2.2: Exemplification of the mechanism of the drift wave (a) and the drift wave insta- 
bility (b). 



potential and the electron density fluctuations are in phase. In the case of perturbations of the 
electric potential, an oscillation of the wave is produced, propagating along the perpendicular 
y-direction with a phase velocity. The drift wave frequency is actually of the order of the 
electron diamagnetic frequency, following the dispersion relation: 

^ ^ ^* ^ _hPlYlI^ (2.24) 

Retaining a possible phase shift between the density and potential fluctuations, turns into 
assuming a response that can be written as Sn^/n = (1 — iSk) e64>k/Te ( k is the wave vector 
in the direction transverse both to the density gradient and to the magnetic field ). This 
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situation is exemplified in Figure [Z2) the radial component of the E x B drift results in a net 
average motion due to the non-zero phase shift between Sn^ and S(f>, in a way that an initial 
perturbation will be sustained and amplified. The same representation illustrates that the 
modes characterized by a positive phase shift, i.e. with Sk < 0, will be damped. In tokamak 
plasmas, there are mainly two possible mechanisms responsible for a non-vanishing phase 
shift: the wave-particle kinetic resonances and the breaking of the hypothesis of electrons 



adiabadicity (cf Eq. (2.21 1). Both of them will be treated in the following of this work. The 



drift-waves represent then an instability mechanism that does not depend on the toroidal 
plasma geometry: for this reason the drift-wave instability defines the so called 'slab' branch 
of several tokamak unstable modes. 

The other big class of instabilities present in the core of tokamak plasmas, are the pressure- 
driven family of interchange modes. The physical mechanism relies on the amount of free 
energy that is released, under certain conditions, by the exchange of two flux tubes around 
a given field line. Formally, this kind of plasma instability is analogous to the hydrodynamic 
Rayleigh-Bcnard instability, whose origin derives from the fact that the fluid temperature 
gradient is aligned with the gravitational force. In tokamak plasmas, the interchange instabil- 
ity is due to both the inhomogeneity of the magnetic field (analogy with the gravity) and on 
the departures from the thermodynamical equilibrium through the presence of large pressure 
gradients (analogy with the temperature gradients). An interchange mode is then unstable 
only when the magnetic curvature, i.e. V-B, is aligned with the pressure gradient Vp. This 
condition is verified only on the low-field side of the toroidal plasma geometry, while the 
high-field side is stable with respect to the interchange drive: the stability condition may in 
fact be written as Vp • < 0. The electric potential is the analogue of the hydrodynamic 
stream function for the Rayleigh-Benard problem. When the condition Vp ■ VB > is ver- 
ified, the E X B, the WB and the curvature drifts combine resulting into a destabilization 
of the convective cells, i.e. iso-contours of the electric potential perturbations. Referring to 





j^^^potential 

(convection cells) 




elecVons 



Figure 2.3: Exemplification of the mechanism of the interchange instability on the tokamak 
low-field side. 



scheme of Figure [273} it appears that the electric drift ve changes of sign between the two 
cells, while the curvature and the VB drifts are vertical, but in opposite directions for the 
ions and the electrons. Globally, the particle motion leads to a positive charge accumulation 
in the positive cell and vice-versa, thus sustaining the electric potential and the convective 
cells instability. 
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The analogy between the tokamak interchange and the hydrodynamic Rayleigh-Benard insta- 
bility is nevertheless limited by a relevant factor. This is the intrinsic anisotropy of tokamak 
micro-turbulence: the unstable modes tend to be aligned to the magnetic field, with a typical 
transverse size remarkably smaller than the parallel one. The tokamak turbulence is then 
often considered as quasi-2D. The fast particle motion along the magnetic field lines make 
that they experience both unstable and stable regions with respect to the interchange drive, 
respectively on low and high field sides: the parallel current appears to be stabilizing in this 
sense. Since the interchange instability can be excited only in presence of a non-vanishing 
magnetic curvature, intrinsic in the tokamak configuration, the related unstable modes define 
a so called 'toroidal' branch. 




Figure 2.4: A typical spectrum of the linear growth rates for the ITG, TEM and ETG 
tokamak plasma instabilities. 



At this point it is useful to introduce the main classes of tokamak instabilities which are 
believed to be the main responsible for the anomalous transport of energy and particles in the 
plasma core. Most of them can be ascribed to the drift-wave or the interchange mechanisms. 
During the motion along the magnetic field lines, the trapped particles undergo the vertical 
drift which is at the origin of the banana shape of their trajectories. The radial extension of 
the banana width Sb can be estimated using the conservation property of the toroidal kinetic 
moment M = mRv^ + eip, giving Sh = -^£=^ where pc is the Larmor gyro-radius. For 

typical values of tokamak safety factor q and local aspect ratio r/R the following ordering is 
obtained: 



(2.25) 



Relation (2.25) can be used to infer the characteristic lengths of the corresponding tokamak 
instabilities, i.e. in the order, TIM (trapped ion modes), ITG (ion temperature gradient), 
TEM (trapped electron modes) and ETG (electron temperature gradient). 



• ITG, Ion Temperature Gradient: these are electrostatic ion modes. They are also often 
referred to as rji modes, since the most relevant parameter for their turbulent drive is 
the ratio r/i = d(logTi)/d(logn). ITG modes include both interchange and drift-wave 
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instabilities, namely (1) slab modes, (2) toroidal modes and (3) trapped ion modes. 
Their characteristic wavelength is bigger or comparable than the ion Larmor radius, 
such that kg Pi < 1.0. 

• TEM, Trapped Electron Mode: these are electrostatic modes due to trapped electrons 
active at ion spatial scales. Their free energy source can be VTg as well Vrig and 
they are usually distinguished between (1) collisionless and (2) dissipative (due to 
collisions) trapped electron modes. The precise limit on their characteristic wavelength 
is somewhat arbitrary, since they largely overlap with ITG modes, but a reasonable 
criterion can be expressed as kgSe ^ 1-0 (where 6^ is the typical radial width of the 
trajectories of the trapped electrons). 

• ETG, Electron Temperature Gradient: these are electrostatic electron modes analogous 
to the ionic rji ones. Their free energy sources are again both VTg and Vrie, i.e. 
Tje — d(logTe)/d(logn), but their typical wavelengths are on the electron rather than 
the ion scales, then kgpe < 1.0. 

• Electromagnetic modes: these can be drift-wave or interchanges instabilities driven by 
the fluctuations of the full electromagnetic potential or micro tearing modes. 

• Fluid-like modes: these can unstable modes active in the plasma periphery and driven 
by Vp, as the resistive balloning modes. 

2.2 Frameworks for the tokamak micro-turbulence 

In the last section, we have introduced the general framework and the basic mechanisms 
which are at the origin of the anomalous transport of energy and particles in the core of 
tokamak plasmas. In order to progress towards a formulation of a transport model able to 
gain reliable predictions, it is important to detail the plasma representation. Historically, the 
plasma dynamics has been described either within a fluid or a kinetic representation. The 
aim of the present section is to justify the relevant reasons for which the quasi-linear model 
here proposed, adopts a kinetic description. The latter is in fact a unique feature among all 
the other actual quasi-linear tokamak transport models, which are instead mainly based on 
a fluid description. 

2.2.1 Kinetic and fluid approaches 

One of the major concern in the representation of the plasma dynamics is the issue of co- 
herence: this kind of constraint is at the origin of relevant increase of complexity of the 
plasma turbulence with respect to more common neutral fluids turbulence. The motion of 
the charged particles in fact, induce electromagnetic fields which back-reacts on the charge 
and current particle densities according to the Maxwell equations. 

A first level of description could try to directly deal with the time evolution of every single 
i particle trajectory, according to the dynamics equation Vi(t) = qs/rUg (E™ +Vi{t) x B™). 
Here the microscopic electric E™ and magnetic fields B™ have to be simultaneously coherent 
to the positions and the velocity of the particles themselves through the Maxwell equations. 
The treatment of this kind of problem, which is at the origin of the Klimontovich equation 
[78], appears as completely untractable since it would end up in tracking typical numbers of 
10^^ particles. 

The starting point of the kinetic representation is the introduction of a phase space which 
contains a significant statistical information of the given system. The six-dimensional phase 
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space takes the form (q,p), where q = (qi,q2,q3) is a three-dimensional spatial basis, and 
p — (pi,P2,P3) is a momentum basis. A distribution function / is then introduced in oder 
to represent the probability density for the system to be found in the elementary phase 
space volume dqdp, i.e. / (q, p, t) = {N (q, p, i)), where the brackets represent an ensemble 
average. Within the Hamiltonian mechanics formulation, it can be shown [78] that the parti- 
cle approach following from the Klimontovich equation, lead to a master Boltzmann kinetic 
equation in the form: 

f - (q • + P • ap) / = |[ - [if, /] ^ c (/) (2.26) 

where the right term C (/) is a generic collision operator, and the general expression for the 
Hamiltonian is: 

H = ^mi| + ecj) + (2.27) 



Depending on the form of this collision operator, Eq. ( 2.26 1 is also referred to as the Fokker- 



Plank equation. Concerning the application of Eq. (2.261 to tokamak plasmas, two relevant 
issues has to be stated: 

• When dealing with hot thermonuclear grade plasmas, the mean free path of the particles 
is usually very large and it can be found that the collision rate is very small. A widely 



used approximation of Eq. (2.26) relies then on considering the coUisionless problem 

the so called Vlasov equation 

[i/,/]=0 (2.28) 



Af /dt — C (/) = 0: this hypothesis defines the so called Vlasov equation 

d/ df 



At dt 



The system defined though the full kinetic equation (2.26) is still of almost untractable 
complexity, since it involves a 6xM dimensional space (where M is the number of 
species in the plasma; often M=3, i.e. ions, electrons and one impurity). In the case of 



strongly magnetized plasma, a further significant simplification of the kinetic (2.26) or 



Vlasov (2.281 equations is the gyrokinetic approximation. Essentially this relies on a 
scale separation argument, averaging the fast gyromotion of the charged particles along 
the field lines, allowing to pass from a6xMtoa5xM dimensional space. 

Even if the kinetic, and more in particular the gyrokinetic approximation, is of great interest 
in magnetic fusion plasmas and will be extensively used in this work, this approach is still 
far to be trivial both from the analytical and also from the numerical point of view. This 
complexity has been at the origin of another major step in the simplification process, defining 
a fluid approach. Here, only the hierarchy of moments of the kinetic equation ( 2.26[ ) is 



considered, through a projection of the latter expression on the velocity basis (l, v, v''...v'') . 
It immediately appears that, thanks to the integration over the velocity space dv^, the 
advantage of this kind of representation is the further dimensional space reduction from 
6xM to 3xM. Consequently, this kind of approach allows a more tractable formulation from 
both the analytical and numerical point of view, allowing to more easily manage the equations 
and to understand the physical mechanisms at play. Historically this is the reason why the 
most part of the actual turbulent transport models deal with a more or less advanced fluid 
description. Nevertheless, two main drawbacks are intrinsic in the fluid approach: 

• Due to the velocity space integration, the fluid moments can hardly account for the 
interaction between waves and particles, as long as resonances in the velocity space are 
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present: the most relevant example is the mechanism of the linear Landau damping. 
This is particularly true for the hot and nearly collisionless thermonuclear plasmas. 
For the same reason, the fluid approach has difficulties to distinguish between passing, 
trapped and suprathermal particles that characterize tokamak plasmas, as well to treat 
the finite Larmor radius effects. 

• The hierarchy of fluid equations obtained by higher order moments on the kinetic equa- 
tion is potentially inflnite; practically, a truncation of this hierarchy at a certain order 
in the v'^ moments represents a closure assumption. The latter is a crucial point for 
any fluid model: a large research activity in tokamak plasmas has been focused on the 
improvement of these closures, aiming to recover the most relevant effect of the kinetic 
approach. Even if this area will not be deepened in this work, the most relevant class 
of advanced fluid closure is the so called Gyro-Landau-fluid set of equations, which is 
widely used for a number of both quasi-linear transport model and also nonlinear nu- 
merical codes. The problem of the fluid closure for nearly collisionless plasmas remains 
nowadays an open issue and a subject of active research. 



2.2.2 The gyro-kinetic approximation 

Thanks to quasi-periodicity of the charged particles trajectories in tokamak plasmas, the 
Hamiltonian H of the system can be described in terms of the angle-actions variables (a, J) 
such that the Hamilton equations are respected: 

a=^=nj (2.29) 



The angle-action variables system a, J is here formally defined. The first pair refers to 
cyclotron gyro-motion: 

ai = flit (2.31) 

Ji = --Ai (2.32) 



The second pair is written as: 



Q!2 = Oai (2.34) 
e J 27r J 



J2 = epassecpT— + f i^rnvw (2.35) 




n2=\^(l> . '^^ I (2.36) 

where Ipass is a constant whose value is 1 in the case of passing particles and otherwise, (pT 
is the toroidal flux normalized to 2tt. Moreover here § = J^^ for the passing particles, while 
for the trapped particles / = 2 where Oq is the angle defining the banana bouncing 
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motion. 

Finally the third pair is: 

as = i^3t (2.37) 

Jg = M ^eil) + mRvgc,^ (2.38) 

^3 = Wa,+^passqW^2 (2.39) 

where gc stands for guiding center, ip for the reference magnetic surface around which the 
guiding center is evolved and A = Vgc,± ■ V(/3 — q (ip) Vgc,j_ • + Jv\\ipd^q. 
This formulation allows to easily recover the integrability of the system and the existence of 
three invariants along the Ji trajectories. According to a perturbative approach, it is then 
possible to rewrite the Hamiltonian H and the distribution function / as: 

H {a, J, t) = Ho (J, t)+J2 ^/^n,^ (J) e^("-°-"*) (2.40) 



/ (a, J, t) = fo (J, t) + y] 5/„,^ (J) e'(n— *) 



(2.41) 



where the time evolution of the equilibrium quantities is supposed much slower than the per- 
turbations one. At a reference equilibrium, the Vlasov problem dfo — [Hq, /q] = will define 
a certain distribution function /o, whose form is a priori fixed also by Coulomb collisions: this 
configures the neoclassical theory problems, which describes the coUisional transport apart 
from any turbulent fluctuation and will not deepened here. Within this framework, a local 
Maxwellian will be adopted: 



/o(J,0 



-Ho{3)/T{t) 



(2.42) 



where n and m are respectively the particle density and mass. 

It is of interest to derive a linearised expression for the response of the perturbation 5f 



to small fluctuations. At the first order, the Vlasov equation (2.28) can be then be rewritten 
as: 



(2.43) 



It is important to note that this results follows from the scale separation on the temporal 
dynamics, i.e. using the quasi-stationary condition for the equilibrium distribution function 



dtfo — [Ho,fo] — 0. Using the Fourier spectral decomposition of Eq. (2.40) and (2.41), the 
following important relation for the linear response of Sf is derived: 



Sfn,u (J) 



n • djHo + iO- 



:Shn,u, (J) 



(2.44) 



where a resonance appears at the frequency oj — n ■ OjHq = n • Hqj- The presence of this 
singularity is at the origin of the Landau damping, a crucial kinetic mechanism which implies 
an energy transfer between the waves and the particles. Formally, the expression does still 
have a validity in the mathematic sense of distributions, when eliminating the indetermination 
on the position of the pole with respect to the real axis. Physically, the solution corresponds 
to a causality constraint, imposing null perturbations for t — > — c»; mathematically, this is 
done through the addition of the imaginary term iO^, where 0^ is small but strictly positive. 
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The presence of the Landau resonance in the hnear kinetic response will be deepened also in 
the following of this work, about the formulation of the quasi-linear transport model. 



A further step with respect to Eq. (2.44 ) can be done within the hypothesis of an equilibrium 



Boltzmann distribution /o of the type of (2.42 1, leading to: 



where 



^/n..(J)=-4^|l- 



w — n • floj + iO 



5hn,oj (J) 



(2.45) 



(2.46) 



Eq. (2.461 is the expression for the generalized diamagnetic frequency, in analogy with the 



definition found by Eq. (2.23). The linear response (2.45) is particularly meaningful from 



the physical point of view: the first term corresponds to the adiabatic response of the dis- 
tribution function, while the second one represents the non-adiabatic component. It will be 
shown that only the latter term originates a turbulent transport via the fiuctuations. 

The gyrokinetic theory was introduced in the attempt of describing the strongly mag- 
netized plasma dynamics over time scales that set apart from the fast gyromotion. In the 
case of tokamak plasmas, where the spatial and temporal variations of the external magnetic 
field are weak, the charged particles undergo 3 types of quasi-periodic motions: (1) the fast 
gyromotion along the magnetic field lines, (2) an intermediate bounce motion along the par- 
allel direction due to the parallel gradients and (3) a slow motion across the field lines driven 
by magnetic curvature and the transverse gradients. Formally, the gyrokinetic theory is an 
improvement with respect to the gyro-center coordinates, introducing a new set of gyrocenter 
coordinates which account for a gyro- averaged perturbed dynamics. A complete overview of 
the modern formulation of the gyrokinetic theory can be found in Ref. [K]. 

The traditional derivation of the gyrokinetic approach follows from a number of funda- 
mental assumptions on the spatial and temporal ordering, comparable to the ones leading to 
the adiabatic theory. The parameter e, ratio between the particle (ion) Larmor radius and 
a macroscopic length L (plasma radius or the density gradient length n/ \V r'n\) is defined: 



L 



< 1 



(2.47) 



Three basic frequency scales exist: the fast cyclotron motion at the frequency lOc, a medium 
frequency which is the typical one for the turbulent fluctuations ujfi ^ ^ euic and the 

slow frequency of the macroscopic transport Wj,. ^ -^e ^ e^ujc- 

The distribution function and the fields are split into a slowly varying (in time and space) 
equilibrium part and a fast varying fluctuating part, such that 



/ (x, V, t) = /o (x, V, t) + Sf (x, V, t) 

The gyrokinetic ordering foresees: 

1. Small fluctuations, of the order of e 

SI \6B\ 
fo ^ |Bo| 



B (x, t) ^ Bo (x, t) + SB (x, t) 
E(x,i) = (5E(x,i) 



\5E\ 



(2.48) 
(2.49) 



\vth'Bo\ 



(2.50) 
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2. Slowly varying equilibrium, on the macroscopic length and time scales (r^r 
the macroscopic transport time scale) 



IS 
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VBo^ 


Bo 
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dt 
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dt 
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3. Fast spatial fluctuations across Bq, on the microscopic length pc 
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4. Slowly varying fluctuations along Bo, on the macroscopic length L 
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(2.51) 
(2.52) 

(2.53) 

(2.54) 



5. Fluctuations varying on a medium time scale (the fluctuations time scale is Tfi 



m^si 

dt Tfi 



dSB 

"df 



SB 



dm 



5E 



6. Collisions act on the medium time scale of the turbulent fluctuations 



'UJfl 



hence 



< 1 



(2.55) 



(2.56) 



These ordering are then defined splitting the distribution function and the fields as a slowly 
varying (in time and space) equilibrium part, characterizing the background quantities, and 
the fast varying fluctuating parts. Since the information on the fast gyromotion dynamics is 
not relevant in this framwork, it is of interest to define an averaging procedure in the Fourier 
space: 



27r 



d0' 

2^' 



,k p. _ 



= Jo {k±Ps) 



(2.57) 



where Jo is the Bessel function. Finally, within the gyrokinetic framework, the general Vlasov 
equation (2.281 takes the following form for a given plasma species s: 

BL _ _ _ 

+ (vexB + Vvb +Vc) • V_l/s +W||V||/3 + /s = 



dt 



(2.58) 



where the hat corresponds to a gyro-averaged quantity. As already mentioned, in order to 
describe a coherent problem, the plasma dynamics is constrained also by the electromag- 
netic Maxwell equations. The Debye length is much smaller than the scale lengths of the 
fluctuations which are described: for this reason, a local electro-neutrality condition can be 
imposed. Within this limit, Eq. (2.58) has to be consistently solved with the Poisson- Ampere 
equations 



(2.59) 
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The solution of the time-evolving nonlinear system described by the set of equations 



( 2.58 )-( 2.59) (mapped on a realistic toroidal magnetic geometry), could in principle provide 
an accurate information about the turbulence dynamics in tokamak plasmas, according to 
the degree of approximation adopted by the gyrokinetic formulation. Nevertheless, the com- 
plexity of the problem is not only analytically largely untractable, but nowadays, and at least 
for the coming 10 years, also too costly from the numerical point of view. 

Solving the linear gyrokinetic dispersion relation 

Solving the linear gyrokinetic plasma dispersion relation, i.e. the linearized gyrokinetic 
equation coupled to the quasi-neutrality condition, is a much easier task with respect to 



the nonlinear system defined by Eqs. (2.58)-(2.59|. Still this approach can provide a great 



amount of information on the tokamak plasma turbulence, hence this is the strategy which 
has been followed by this work. 

This paragraph is dedicated firstly to illustrate the general structure of the gyrokinetic elec- 
trostatic linear problem, and therefore a particular formulation that is used to compute 
numerical solutions. Moreover, the following discussion will provide the basic framework 
which is employed by the quasi-linear transport model QuaLiKiz. QuaLiKiz in fact is en- 
tirely based on a linear gyrokinetic eigenvalue code, Kinezero jl2) . The hypotheses and the 
approximations underlying Kinezero are then completely shared with QuaLiKiz. 



The electro-neutrality constraint appearing in the first of Eqs. (2.59) can be rewritten 
according to a variational approach [38! , giving: 

(uj) = with C,{Lo) = -Y,j eMl,,^ (x) <5</.;,^ (x) (2.60) 

s n.u) 

The linear gyrokinetic plasma dispersion relation is obtained combining the linearized re- 



sponse for the distribution function Eq. (2.45), with the quasi-neutrality condition Eq. 



(2.60), leading to: 



{ 



n VI* 



n • n 



OJ 



^0^ 



a, J 



= (2.61) 



Formally, the solution of the Eq. (2.61 ) provides for each n one or more complex eigenvalues 
u) = Ur + i^. If 7 > 0, then the eigenvalue is associated with a linear unstable eigenmode 
which is exponentially growing in time with the finite real frequency Ur- 



The solution of the eigenvalue problem defined by Eq. (2.61) is still not trivial even from 



the numerical point of view, because of the high order of the matrix which are involved. Re- 



cently, accurate eigenvalue numerical solvers for the Eq. (2.61 ) have been developed [571 155] : 
even if less expensive than the full nonlinear simulations, presently these codes still require a 
significant amount of computational resources, that make them not immediately applicable 
within a time-evolving integrated transport solver. Consequently, the quasi-linear transport 



model here proposed is based on an further approximation of the Eq. (2.61), which is nu- 



merically solved by the code Kinezero. One of the most relevant feature of this code is in 
fact its significantly lower computational requirement with respect to the other gyrokinetic 
linear solvers. The approximations that are employed by Kinezero are then here below briefly 
recalled. 



1. Lowest order of the ballooning representation. In the coordinate system (r, 9, (p), the hy- 
pothesis of perfect toroidal axisymmetry implies that e™"^ is an eigenvector. Moreover, 
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the tokamak micro-instabilities are characterized by the strong anisotropy such that 
fc|| k±. The ballooning representation uses these arguments in order to reduce the 
spatial dimensionality of the problem [53]. In particular, the eigenvector S4> {r^0,(p,t) 
is re-written as: 

6cl> (r, 0.^,t)^Y. j^H ^'^-^o {0 -f 2^1) e^M^-<i(r)(e+2.i-e,))-.t] (2.62) 

n.u) 1— — OQ 

The ballooning angle 6*0 is here taken to be 0^ — 0, since in the most cases the unstable 
modes are ballooned on the low-field side. Within the lowest order of the ballooning 
representation, the second radial derivatives of the equilibrium quantities are neglected, 
i.e. \d\ < (VrlogA) ^, where A stands for a generic equilibrium quantity and d is the 
distance between two resonant surfaces d = — (nVr^) . This lowest order limit is 
formally no longer valid when the first radial derivative of an equilibrium quantity is 
approaching to zero: this applies in particular when the magnetic shear \s\ — 0. 



2. Trial Gaussian eigenfanction. In order to quickly find the eigenvalues of the linear 
dispersion relation, Kinezero adopts a trial S(f> eigenfunction. The latter one is chosen 
to be the most unstable analytical solution in the fluid limit for strongly ballooned 
modes. It can be shown that S(f> (kr) is a Gaussian in the form: 

<50(fc,) = J0o^e-^>'/2 (2.63) 

The mode width w is an important parameter which is as well calculated in the fluid 
limit, assuming the interchange as the dominant instability. The complete expression 
for w can be found in Ref. |12j . 



3. Functionals for trapped and passing particles. The complex zeros of the Eq. ( 2.61 1, re- 
written as V (w) = J2s "^5^ ~ (^)tr s ~ i^)ps sj ~ ^ ^^'^ numerically solved, where 
the ion, electron and one impurity species are accounted in both their trapped tr and 
passing ps domains. The frequencies are: 



keTs 
e,BR 



£- 



Js 



keTs 
e,BR 



£ (2 - Xb) f (X) + knvn 



(2.64) 
(2.65) 



where kg — nq (r) /r is the toroidal wave-number and / (A) a function of A depending 
on the magnetic geometry and the MHD parameter a, differing for trapped and passing 
particles jH] . In the resonance frequency both the vertical drift and the parallel motion 
appear, while the terms eventually coming from V?;|| and Vw_l are neglected. 
Firstly focusing on the trapped particles, the corresponding velocity integration domain, 
using the £, A variables, is: 



/£,A 



trapped 
1 



+ 00 



d£ I dX-^jQ {k±pcs) Jq {krSs) 
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^i-x(T,e)b 



HB (r, 9 = 0) 
fiB {r, e^n) 



(2.66) 
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Jq (fcj_pcs) and Jq {krSs) are Bessel functions standing respectively for cyclotron gy- 
roaverage and the bounce average. As the trapped particles are mainly efficient at low 
wave-numbers, the hypothesis that n <^ pv^j therefore krd tt is used. Moreover, 
because of the assumption of a local MaxweUian equilibrium, the energy integration of 
the Bessel functions are separated from the rest. Finally, the collisions on the trapped 
electrons are also retained through a Krook operator |57j . 
Concerning the passing particles velocity integration, this can be written as: 



/ -^d£ dX—-Y,^jSik±Pcs) (2.67) 
Jo Jo 4a;2 2 ^ 



where the hypothesis of n >• p^^, therefore krd <C tt, is used. Also the A integration 
is approximated, so that: 

nfls ~ (nils);^ ~ (^IPIi)a ~ — -j^wVS (2.68) 

A wide documentation on the benchmark efforts made comparing Kinezero against 
other linear gyrokinetic codes is reported in Refs. ^12^ S2i ilSj . 



2.2.3 Choosing the kinetic approach 

In the previous paragraphs we have briefly introduced the two most common frameworks 
used in the description of the tokamak plasma turbulence, namely the kinetic and the fluid 
approach. A unique feature of the quasi-linear model proposed in this work is the use of 
a gyrokinetic formulation, while most of the actual quasi-linear transport model are fluid 
|104| or gyro-fluid |103[ 157] . This paragraph and the next one will be devoted to provide 
demonstrations that supports the preference in favor of the kinetic approach. In particular 
it will be shown that: 

• The fluid approach, which intrinsically requires an arbitrary closure, typically overes- 
timates the values of the linear threshold of the tokamak micro-instabilities. 

• The kinetic framework consistently recovers the fluid limit. 

• Several mechanisms which are at play in tokamak plasmas are a direct consequence of 
the presence of the kinetic resonances. 

A comprehensive understanding of the fluid approximation is based on the derivation 
following from the kinetic approach. The main interest of a kinetic treatment for deriving 
the linear unstable modes thresholds concerns the accurate treatment of the wave-particle 
resonance at the frequency a;^ such that: 

= djHo ~ LOD (wf,wi) + k\\v\\ (2.69) 

where ujjj is the magnetic (curvature and VB) drift frequency, fcy is the parallel wave vector, 
while W|| and v± are, respectively, the parallel and perpendicular velocities. We can assume 
that the magnetic drift component of the resonance will be dominant on the instability growth 
rates, since the parallel dynamics is subdominant for s/g ^ 2rl Neglecting the term in 
the resonance and the finite Larmor radius (FLR) effects, the linear kinetic density response 



'^s is the so called magnetic shear, defined as s = rVrg/<? 
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for the s species Sug directly follows from the linear response (2.45) and can be written in 
the following way: 



(2.70) 



where the characteristic frequencies are expressed as lu* — w*^ + uj^^ {£ — 3/2), lujjs = 
\sUJdts£ and As = cosx + sxsinx (% is a coordinate along the magnetic field line). The 
resonance and the diamagnetic frequencies are analogous to the ones already introduced by 



Eqs. (2.64|-(2.64|; here the vertical drift and the diamagnetic terms driven by the density 



and the temperature gradients are distinguished, such that: 



WDTs — 2fc6 



e,B B 



(2.71) 



' CsB ns 



Ts VrTs 
BsB Ts 



(2.72) 



where kg is the poloidal wave-number. 

Dealing firstly with ion ITG modes, a very crude fluid limit (without any closure) can 
be derived from the linear kinetic density response through a second order expansion of the 
Eq. (2.70) based on the condition ujoi/^ *C 1, i-e- considering that the mode frequency 



is far from the resonance. The resulting ion response Sni/ui has then to be coupled to the 
quasi-neutrallity condition Sni — Sn^ and assuming adiabatic electrons {5n(,/ne = ed(f)/Te). 
The following second order dispersion relation is then derived: 
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(2.73) 



Since at the threshold the imaginary part of oj is vanishing, the following expression for the 
critical ITG R/Lxi can be obtained: 
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(2.74) 



It is interesting to compare the latter expression for the ITG threshold obtained by 
a primitive fluid expansion, with a more sophisticated fluid model which accounts for a 
particular closure. Leaving the derivation details in the appendix [Aj the advanced fluid 
Weiland model ^104j operates through this strategy, providing the following expression for 
the ITG linear threshold: 
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(2.75) 



When comparing Eq. (2.74 1 and Eq. (2.751, it appears that a relevant difference is carried 
by the parametric dependence of the thresholds on the ratio Ti/T,: in particular, while the 
Weiland ITG threshold increases for higher Ti/T,, the opposite behavior is shown by the 
crude fluid limit. In fact, as it will be demonstrated in the following, for conditions close 
to the flat density limit, the hypothesis of mode frequencies far from the kinetic resonance 
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(crude fluid limit) completely fails in reproducing the correct Tj/Te dependence of the ITG 
threshold, which is instead linked to kinetic resonance effects. The closure introduced by 
the advanced fluid models has specifically the aim of recovering some of the essential physics 
carried by the presence of the kinetic resonance. The linear threshold derived from the 
Weiland advanced fluid model appears then to be a good candidate to be compared with 
completely kinetic approach. Before this direct comparison, a brief discussion on the linear 
threshold of the trapped electron modes is here given. 

The problem of the VTg linear instability threshold for trapped electron modes has been 
quite extensively addressed through both analytical efforts and numerical simulations. The 
TEM threshold problem can be again addressed starting from the linear kinetic response 
given by Eq. (2.70) in the presence of a nonzero fraction of trapped electrons ft In this 
case, the hypothesis of retaining only the magnetic drift contribution of the kinetic resonance 
is even more justified, since passing electrons are considered adiabatic. On the other hand, 
a more subtle question is arising when coupling the electron kinetic response with the ion 
one. In fact, if for ITG modes the hypothesis of electron adiabaticity can be valid at high 
collisionality, the symmetric choice of adiabatic ions is not feasible, since i^ei 3> (where 
Vij is the collision frequency between species i and j). Hence, the non-adiabatic response of 
both trapped electrons and ions through quasi-neutrality is retained, giving: 



I -ft 



Te 



1 - 



(2.76) 



The choice for the ion response is put in evidence through the parameter H appearing in the 



last term of Eq. (2.76). Even if H will depend primarily on the modes frequency w and on 



both and VrT, two extreme cases can be drawn: keeping the condition H = 1 means 
adopting ion adiabaticity, while on the other hand H = implies entirely neglecting the ion 
response. 



Even in this case, it is initially possible to adopt a crude fluid limit based on Eq. (2.76), 



considering conditions far from the electron resonance, i.e. ujuel^ <C 1; a second order expan- 
sion for the non-adiabatic electron response is then allowed. According to this approximation, 
the following R/Lxe TEM threshold expression can be derived: 
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where Kt 



and T - 



Within this simple fluid approximation, the two opposite 



ft 

limits H = 1 and H = lead to discrepancies in the temperature ratio behavior of the TEM 
threshold. In fact, when choosing H = 1 (adiabatic ions), a TEM R/Lxe threshold raising 
with higher Ti/T^ is found. Within the second limit H = instead, the TEM threshold 
exhibits no temperature ratio dependence, as often found in literature |104|, I83j. In fact, as 
done for the ITG threshold, Eq. (2.77) can be compared with the prediction of the Weiland 
advanced fluid model that gives: 
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(2.78) 



Finally, it is relevant to quantitatively compare the expectations for the linear ITG and 
TEM thresholds from the kinetic and the advanced fluid approaches. In Fig. |2.5[ the para- 
metric space referred to the normalized density and temperature gradients (R/ Ln, R/ Lt) is 



^In this case the quasi-neutrahty condition will be modified as Sui 



■ 5n 



e,trap 



where the non- 



adiabatic trapped electron response is included, while passing electrons are still treated as adiabatic 



24 



Description of the plasma instabilities 




Figure 2.5: (a) ITG threshold according to the Weiland advanced fluid model (full lines) and 
computed by the linear gyrokinetic code Kinezero (dotted lines) in the plane R/ LmR/ Lt 
for different values of the ratio Ti/T^- (b) TEM threshold according to the Weiland advanced 
fluid model (full line) and computed by the linear gyrokinetic code Kinezero (dotted lines) 
in the plane R/Ln, R/Lt for different values of the ratio Ti/T^. 



explored, distinguishing between linearly stable and unstable regions according to the Wei- 
land model and to the numerical results obtained using the linear gyrokinetic code Kinezero. 
It clearly appears that even the advanced fluid approach overestimates the linear threshold 
with respect to the fully kinetic result. Even if the closure is able of reproducing the overall 
Ti/Tf. dependence of the ITG threshold, the critical gradients significantly differ between the 
two formulations. Moreover, no temperatures ratio dependence for the TEM threshold is 
expected according to the Weiland fluid model, while the linear gyrokinetic simulations re- 
veal that also the TEM threshold is affected by the Ti/T^ variations. The fluid overestimates 
would then imply non-negligible errors when applied into a quasi-linear transport model, 
since no turbulent transport can be predicted in the linearly stable region. Therefore, there 
are strong indications to prefer the kinetic framework if one hopes to realistically describe 
the tokamak turbulent transport using a quasi-linear theory. 



2.2.4 The T^/Te dependence of the hnear ITG-TEM thresholds 

In order to reinforce the argument in favor of the use of the kinetic approach with respect to 
the fluid one, a particular issue which is relevant for realistic tokamak applications, will be 
here below discussed in detail. As anticipated in the previous paragraph, this example deals 
with the temperatures ratio T^/Te dependence for the linear threshold of the ITG-TEM un- 
stable modes. Firstly, a simple analytical derivation will be used to highlight that the correct 
temperatures ratio dependence of both the ITG and TEM linear thresholds is directly linked 
to the presence of the kinetic resonance, providing expectations in agreement with the lin- 
ear gyrokinetic simulations. Secondly, it will be shown that the fluid expectations naturally 
appear as a particular limit of the kinetic framework. The whole analysis help moreover in 
understanding the impact of the relevant parameter T^/Tg on the tokamak micro-turbulence. 

From the experimental point of view, a large number of evidences have revealed that 
the plasma performance is significantly improved when the ratio between ion and electron 
temperature Ti/Tf. is increased [1061 ISl [1091 [711 [H^ lllj . The hot ion high confinement mode 
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(H-mode) has allowed achieving the highest fusion yields in the Joint European Torus (JET) 
IHHinS]; similar results have been reported also in DIII-D [T7], with the 'supershot' in the 
Tokamak Fusion Test Reactor (TFTR) j53j and with the advanced scenarios in ASDEX- 
Upgrade [90]. However in the next step devices, the dominant electron heating by fusion 
a-particles in the center coupled with the thermal equilibration between ions and electrons 
at high densities, will lead to a ratio Ti/Tg. slightly below unity. Hence, it is of interest 
to have a theoretically based understanding about the impact of the temperatures ratio on 
micro-turbulence, which strongly affects the plasma performance. 

Some dedicated transport analysis have already assessed the Ti/T^ dependence on the 
energy confinement time O |84] and on the ion heat transport together with the role of 
the radial gradient of toroidal velocity Vvtor [12] ■ Nevertheless few systematic analytical 
or numerical study of the temperatures ratio impact on the instability thresholds of ion 
temperature gradient ITG modes and trapped electron modes TEM has been carried out. 
Experimentally, the existence of such a threshold on the electron temperature gradient length 
has been proved [52( I88j . Several formulations for the critical normalized temperature gra- 
dient length R/Lti.c — R\'^rTi^e\/Ti^e of these instabilities have been proposed according to 
the fluid or the kinetic approaches. In the case of pure ion modes, which have been firstly 
thoroughly studied, a threshold increasing with Ti/Te is well supported inside all the present 
formulations [451 1104j . Conversely, no impact coming from temperatures ratio variations was 
foreseen for the electron modes threshold |104l IBS] . This study has systematically assessed 
for the first time also the Ti /Te dependence of the TEM threshold [21] . 



TEM threshold within a kinetic approach 



Firstly, for reasons of simplicity, pure trapped electron modes in the strict flat density limit 
will be considered, thus implying R/Lxi — R/L^ — 0; these assumptions will be later relaxed. 
The problem is addressed starting from the kinetic dispersion relation (2.76) considering 
a mode in the electron diamagnetic direction. The threshold for the mode instability is 
derived isolating the imaginary contribution coming from the resonance, using the relation 
hm£_>o ^ = PP (s) T iTi-^ (x), thus giving: 
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At the threshold the imaginary part of Eq. (2.79) is vanishing, leading to 



''^e 7, 



(2.80) 



When using the latter relation into the real part of Eq. (2.79) (note that the dependence on 
the normalized energy £ will not be present anymore) , the following relation can be obtained: 
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(2.81) 



A last approximation regards neglecting the integration over A. Even if A accounts for the 
mode structure through a x coordinate along the field lines, if the eigenfunction has a strongly 
ballooned structure, lodtc can be replaced with its value at x « [HHl US]. Within these 
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hypotheses, Eq. (2.81) reduces to a basic expression for the TEM instability threshold: 
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(2.82) 



where the dependence on the parameter Ae has been left explicit. In the case of trapped 
electron modes, the localization of the eigenfunction can be affected by the minor role of 
parallel dynamics; for this reason the value of Ae is not expected equal to 1 and it will be 
treated as a free parameter including a dependence on the magnetic shear s |36| . 



Equation (2.82) highlights that explicit T^/Ti dependence for the R/Ltc TEM threshold 
can be potentially achieved for H ^ 0. An implicit expression for H has been analytically 



obtained and numerically solved making use of the resonant condition (2.801: the details on 



this derivation are discussed in Appendix A. Equation ( 2.82 ) then provides us a new reference 
for studying the R/Lxe TEM threshold scaling with T^/Ti in the flat density limit. 

This analytical approach has been tested against linear gyrokinetic simulations with the 
code Kinezero. In the framework of the Ti/Tg study here presented, the set of plasma pa- 
rameters reported in the Table [2?T] is defined in order to test the analytical predictions. 





Ro/a 


r/a 


q 


s 


CtMHD 






2.8 


3.0 


0.4 


1.4 


0.8 


0.0 


1.0 


0.0 



Table 2.1: Plasma parameters adopted for the linear gyrokinetic simulations here presented 
and performed with the code Kinezero. 

The numerical simulations consider 30 toroidal wave numbers in the range 0.1 < kgps < 2.0. 
CoUisionless trapped electron modes in the flat density limit have been simulated imposing 
R/Lxi = and R/Ln = 0; a scan on the ratio Te/Ti from 0.25 to 4 has been performed fixing 
alternately T^^i = 4 keV. The instability threshold has been identified within the interval of 



R/ Lxe values where the linear growth rates become nonzero. Results are plotted in Fig. 2.6 
together with the analytical predictions for the TEM threshold and for H. With the present 
set of plasma parameters the value of Ag = 0.66 has been considered; the latter choice is 
reasonable if compared to the form proposed in Ref. [35] giving Ag = 1/4 + 2s/3 = 0.783. 

Several new and interesting features emerge from results shown in Fig. |2.6| First of all, 
an effective dependence of the R/Lxe TEM threshold on the electron to ion temperature ratio 
is both numerically and analytically recognized in the flat density limit. This evidence has 
been observed here for the first time and it is particularly remarkable for tokamak relevant 
conditions. Indeed a linear fit for 0.5 < T^/Ti < 1.5 gives in fact R/LreM = 1.1^+4.4. The 
Te/Ti dependence weakens for Tf, ^ Ti where the ion response becomes negligible. Moreover 
this feature is clearly not dependent on fixing Te or Ti and numerical results confirm that 
having neglected FLR effects does not turn into a too severe simplification. 
In the flat density limit, the TEM threshold is therefore largely dominated by the kinetic 
resonant effects. Moreover, consistent evaluation of the non-adiabatic ion response is crucial 
for successfully reproducing the linear increase of TEM R/Lrcth with T^/Ti at V^n — 0, 
especially for T^ ^ Ti. A final statement has to be specified concerning the range of va- 



lidity of the present approach. Equation (2.82) has been in fact derived retaining only the 



resonant contribution of the electron response; this means that the validity of this choice is 
limited to conditions of temperature and density gradients satisfying R/Ln <~ 3/2R/ Lt^ih 
and R/LT,th >~ 2 (see Appendix [A| for details). The problem of the temperatures ratio 
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Figure 2.6: (a) Tg/T,; dependence of the TEM instability threshold calculated with 
KINEZERO in the flat density limit at fixed Ti and Tg. The numerical thresholds present 
intrinsic error bars, whose amplitude simply depends on the chosen step size for the R/Lxi 
scan. Analytical predictions following from Eq. (2.82 ) are also plotted considering Ae = 0.66. 
(b) Analytical estimates for the H parameter appearing in Eq. (2.82). 



threshold dependence beyond this limit requires then a different approach. 

Conditions of both V^ti = and V^Ti — adopted until here are not commonly compat- 
ible with realistic plasma scenarios. In the following, we will discuss the temperature ratio 
dependence of the TEM R/LTe,th in the presence of non-adiabatic ion response driven by 
nonzero V^ri and S/rTi. 

Modes rotating in the electron diamagnetic direction obeying the linear kinetic dispersion 



relation (2.791 are considered. Nonzero R/Lxi will be simply treated as a constant, playing 
the role of driving additional contribution in the non-adiabatic ion response. In other words, 
the possible contemporary presence of two unstable solutions (electron and ion branches) is 
not retained. Within our analytical approach based on the electron resonance condition, the 
ion non-resonant response has been evaluated including the additional contributions of both 



a;*j,w^j 7^ 0, as detailed in Appendix A. In this case the value of H appearing in Eq. (2.791 
depends on the parameters Ti/T^, R/Ln, R/Lxi, and ft- The effect of magnetic shear is 
incorporated in the free parameter Ag. 

This approach has been tested against linear gyrokinetic simulations with Kinezero, per- 
formed on a single wave number corresponding to the maximum of the linear TEM spec- 
trum {kgps ~ 0.5) and looking for the critical normalized gradients. Conditions of domi- 
nant TEM turbulence have been achieved in the simulations imposing R/Lxi = 0.65R/ Lxe^ 
while the analytical calculations have been performed fixing R/Lti = 2.5; finally the values 
Ti/Tf, = 0.5, 1,2 have been considered. Both simulation results and analytical expectations 
for the TEM threshold are presented in Fig. |2.7[ 

A good agreement between analytical and numerical results is observed; small discrep- 
ancies have to be ascribed to the role of R/Lxi, which is kept constant within the analytical 
approach while it is multiplied by a factor of 0.65 with respect to R/LTe,th in the simula- 
tions. At zero density gradients, but in presence of non-negligible VrTi, numerical results and 
analytical expectations agree in recognizing an increase of the R/L^e TEM threshold with 



higher Tg/Ti. This behavior has been already found with R/Ltz = in Fig. 2.6 meaning 
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Figure 2.7: (a) Analytical predictions for the R/Ltc TEM threshold vs R/Lm imposing 
R/Lti = 2.5 for different ratios Ti/T^. (b) R/Ltc TEM instability threshold calculated with 
Kinezero on a single wave number (maximum of the linear spectrum) considering R/Lti — 
0.65R/ Lxe and for the same ratios Ti/T^. 



that a modest ion temperature gradient is not sufficient to significantly affect the T^/Ti TEM 
threshold scaling. 

For higher values of density gradients, above R/L„ ~ 1.2 for this set of plasma parameters, 
an inversion of the temperature ratio dependence of R/LTe,th is observed. In this case, the 
VrTe TEM threshold is slightly increasing with higher Ti/T^. From the analytical point 
of view, the peculiar reversal in the temperature ratio dependence can be explained only if 
accounting for the ion response dynamics within the kinetic modes dispersion relation. This 
effect is clearly due to the role of nonzero V^n in the non-adiabatic ion response. 



ITG threshold within a kinetic approach 



We have anticipated that the correct temperature ratio dependence of the ITG threshold 
in the flat density limit is due to the role of the kinetic resonance. This result was already 
known by previous works [861 and it can be also recovered by advanced fluid model like the 
Weiland model |104| ; in Appendix A the derivation of such a model is briefly recalled, high- 
lighting the presence of a resonant denominator acting within the fluid dispersion relation. 
Here instead, the Ti/T^. dependence of the ITG threshold in the absence of density gradients is 



derived directly from the ion linear kinetic equation (2.70), coupled with the quasi-neutrality 
condition and adiabatic electrons. Similarly to the procedure adopted for the electron modes, 
isolating the contribution coming from the ion kinetic resonance (imaginary part) provides a 
condition on the mode frequency: 

w = <:+Wt. f— -1) (2.83) 



\W£)Ti 2^ 

The latter relation can be used for rewriting the real part of the linear ion kinetic response 



(2.70) coupled to adiabatic electrons, resulting in the following threshold expression: 

R 

Ti th 



L 



1 + ^) (2.84) 



This simple relation has been tested against numerical simulations using KINEZERO with 
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R/Lxe = and R/L^ = 0. A scan on the ratio TijTf, from 0.25 to 3.0 has been performed 
alternately fixing Te.i — 4 keV. The resulting ITG threshold versus Ti/Tf, is presented in Fig. 
2.8 Numerical simulations highlight that the ITG instability threshold is clearly raising with 
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Figure 2.8: Ti/T^ dependence of the ITG instability threshold calculated with KINEZERO 



in the flat density limit at fixed Ti and T^. Analytical predictions following from Eq. (2.841 
are also plotted. 



increasing Ti/T^ on the whole interval. Despite having neglected the parallel ion dynamics 



and simplified the A-integration, the agreement with the analytical form by Eq. (2.84) is 
quite satisfactory. Some discrepancies are nevertheless observed when fixing the ion rather 
than the electron temperature; these difl'erences are most probably ascribed to the role of 
passing ions. In the flat density limit, the ITG instability threshold is largely dominated by 
the role of kinetic resonance (magnetic drift contribution) , which is responsible for its clear 
increase with higher TijT^. 

Electrons are usually retained adiabatically when deriving the ITG threshold; neverthe- 
less including the non-adiabatic response due to trapped electrons (hence nonzero V^n and 
VrTe) leads to a more accurate treatment. As detailed in Appendix A, an analytical proce- 
dure very close to the one already applied for TEM has been adopted. Analytical predictions 
have been tested against numerical results. The ion modes thresholds have been calculated 
with Kinezero on a single wave number (maximum of the linear ITG spectrum) for ion tur- 
bulence dominated plasma. In the simulations, R/Lxe = OAR/Lxi has been imposed, while 



R/Ltc = 3 has been considered within the analytical calculations. Fig. 2.9 summarizes the 
results for TjTe = 0.5, 1, 2. 

The ion modes threshold does not exhibit any inversion in the temperature ratio de- 
pendence, which conserves the same Ti/Tf, scaling found in the flat density limit. Including 
the electron non-resonant response leads to a slight decrease of the ITG threshold at higher 
density gradients, as confirmed by the simulations; the latter effect is not captured by the 
advanced fiuid models. 

Briefiy, kinetic resonant effects playing in the temperature ratio ITG threshold dependence 
have been verified beyond the usual hypothesis of adiabatic electrons. The increase of ITG 
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Figure 2.9: (a) Analytical predictions for the ITG threshold vs R/Ln imposing R/Lte = 3 
for different values of T,/Te. (b) ITG instability threshold calculated with KINEZERO on a 
single wave number (maximum of the linear spectrum) and considering R/L^e = OAR/LTi 
for the same ratios Ti/T^. 



threshold with higher Ti/Tf. is not significantly affected by nonzero V^fi and V rTe- If in the 
strict flat density limit, ion and electron modes thresholds have opposite temperature ratio 
dependence, a critical value in R/Ln exists, above which an analogous scaling of both the 
critical R/Lxi.e increasing with higher Ti/T^ is expected. 

Linear stability diagrams and Ti/T^ impact on ITG-TEM grovifth rates 



The mode frequency relations (respectively Eq. (2.83) for ion and Eq. (2.80) for elec- 
tron modes) are consistent with the hypothesis of modes propagating in the ion (electron 
respectively) diamagnetic direction; the latter condition turns into the following limits: 
R/Ln <« 3/2R/LT,th and R/LT,th >~ 2 (see Appendix A). These criteria define the region 
where resonant kinetic effects are dominant; until here, the X-JT^ dependence of ITG and 
TEM temperature gradient thresholds have been studied within these limits. Here instead we 
extend the study of temperature ratio dependence of instability thresholds when considering 
frequencies away from the magnetic drift ones, where the non-resonant effect will impact 
the thresholds: in other words, the latter conditions correspond to what has been previously 
referred as the crude fluid limit. 

Stability diagrams for ITG and TE modes in the plane {R/ Ln, R/ Lx) have been nu- 
merically obtained with the code Kinezero. Simulations have been performed considering 



R/Lxi = 0.6R/LTe (Fig- 2.10) on a single wave number corresponding to the maximum of 
the linear ITG-TEM spectrum. Two unstable branches coexist, corresponding to modes in 
the ion and in the electron diamagnetic directions. Concerning the ion branch, the resonant 
kinetic effects are dominant at low density gradients and the ITG threshold clearly raises 
with the ratio Ti/T^. At higher Wrn/n and WrTi/Ti, non-resonant terms (and therefore the 
crude fluid limit) become more and more relevant; this results in a reversed Tp/Ti scaling of 



the ITG threshold, as already predicted by the analytical fluid limit of Eq. (2.74). 

For a wide region of the plane {R/Ln, R/Lt), resonant kinetic effects are indeed at play 
for coupled ITG-TEM modes, signiflcantly affecting their instability thresholds. In this case, 
for R/Ln ^ 1.2 these results highlight that raising the ratio Ti/T^ with R/Lti = R/Ltc 
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Figure 2.10: Stability diagrams of both ion and electron modes at different ratios Ti/Tg 
calculated with KINEZERO on a single wave number (maximum of the linear spectrum=, 
(a) imposing R/Lxi = 0.6R/LTe, (b) imposing R/Ln = R/Lxe- 



is simultaneously widening the stability region of both ion and electron modes through a 
direct increase of their instability thresholds. This observation could eventually be regarded 
as a theoretical insight for the large amount of experimental evidence that recognizes the 
beneficial effect of high Ti/Tf. on the plasma confinement. Nevertheless, one of the main 
points raised by this work highlights that the Ti/T^ thresholds scalings are far from being 
universally valid, since they can be reversed at higher values of density gradients. Moreover, 
the results presented in this work are referred to a specific set of plasma parameters (namely 
safety factor g, magnetic shear s, fraction of trapped particles, etc.) and the mode thresholds 
are sensible to all of them; a lot of care is then required when comparing these predictions 
with the experiments. More realistic self-consistency between different transport channels 
(ion/electron particle and heat transport) has not been addressed here and it is beyond the 
scope of this present analysis. 



A final relevant question could be: do these revised temperature ratio thresholds depen- 
dences actually turn into appreciable variations of the modes linear growth rates? 
The growth rates are in fact expected to strongly affect turbulent heat and particle fluxes. 
For this reason the transition from TEM to ITG turbulence has been studied with Kinezero 
using the set of plasma parameters of Table [2TT] The electron temperature gradient has been 
fixed sufficiently high to have TEM turbulence active, while a scan over R/Lxi has allowed 
moving from electron to ion turbulence at Ti/T^ = 0.5, 1, 2. Moreover the role of normalized 
density gradients has been taken into account considering two different cases: R/Ln = 



(Fig. 2.11) and i?/L„ = 3 (Fig. 2.12 1. In these conditions, ITG and TE modes are often 



coupled; the following plots distinguish the growth rates of modes propagating in the ion 
(7oi) and in the electron (7oe) diamagnetic direction. 

Results of Fig. |2.11| in the flat density limit exhibit consistent behavior with both the 
previously found analytical and numerical thresholds scalings. The net effect of raising the 
ratio Ti/T^ has opposite effects on the two branches, i.e., stabilizing for ion modes but desta- 
bilizing for electron ones; moreover the transition from electron to ion modes shifts towards 
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Figure 2.11: Linear growth rates vs R/Lti considering R/Ln = and R/Lxe ~ 6; both the 
first and the second most unstable solution (in the ion and electron diniagnetic directions) 
are plotted for (a) TJT^ = 0.5, (b) Ti/T^ = 1.0, (c) TJT^ ^ 2.0. 




R/L^. R/L^. R/L^. 



Figure 2.12: Linear growth rates vs R/Lti considering R/Ln = 3 and R/Lxe = 5; both the 
first and the second most unstable solution (in the ion and electron dimagnetic directions) 
are plotted for (a) T./T^ = 0.5, (b) Ti/T^ = 1.0, (c) T./T^ = 2.0. 
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higher R/Lxi when increasing Ti/Tg. Hence when moving from TEM to ITG dominated 
plasmas at neghgible Vrn/n, one expects a reversal of the T^/Tg impact on the confinement. 



In the case of Figs. 2.12 the role of nonzero density gradients R/Ln = 3 reverses the temper- 
ature ratio dependence of electron linear growth rates, while leaving unaffected the impact on 
the ion ones. Here in fact the increase of Ti/T^ lowers not only the ion but also the electron 
growth rates, breaking the opposite Ti/T^ scaling observed in the flat density limit. As the 
thresholds, the linear growth rate behavior of both ITG and TEM modes due to the role of the 
temperature ratio can undergo deep changes depending on the values of the density gradients. 

In summary 



1. An analytical approach based on the kinetic resonance contribution for deriving esti- 
mates on the instability thresholds is able of successfully reproduce the results from 
linear gyrokinetic simulations. When considering modes in the electron (ion resp.) 
magnetic drift direction, a relevant point is also the self-consistent treatment of the ion 
(electron resp.) response in both its adiabatic and non-resonant contributions. 

2. At low values of normalized density gradients, the ITG R/LTi,th increasing with higher 
Ti /Te is due to the drift magnetic resonant contribution acting within the kinetic mode 
dispersion relation. Regarding pure trapped electron modes at Vru/n = 0, significant 
raising of R/LTe,th with higher Te/T^ has been analytically and numerically found for 
Ti K, Te- Hence linearly, at low V^n/n the plasma confinement with higher Ti/Tg is 
expected to degrade for TEM turbulence and to improve for ITG turbulence. 

3. When instead increasing the normalized density gradient at VrTi/T,; ^ 0, the T^/Ti 
scaling of the VrTe/Tg TEM threshold reverses, due to the non-adiabatic ion response: 
in these conditions, increasing the ratio Ti/Tg turns in raising both TEM and ITG 
thresholds. Similar inversion is not found for the ion modes threshold. The numerical 
linear gyrokinetic stability diagrams in the plane {R/ L^, R/ L^) allow to identify the 
conditions where resonant effects are overcome by non-resonant terms. 

4. A major step is still required in order to apply the present understanding to the tokamak 
discharges. The self-consistency between particle and heat transport channels for both 
species and the sensibility to several plasma parameters have to be accurately accounted 
for within a physically comprehensive gyrokinetic quasi-linear transport model. 

5. A significant region for feasible tokamak plasmas scenarios in the plane {R/Ln, R/Lt) 
has then been studied and identified as dominated by the physics of magnetic drift 
kinetic resonance. A comprehensive kinetic approach is essential to correctly capture 
both the resonant and the non- resonant components of the plasma response, without the 
need of any additional arbitrary closure. Finally, the gyro-kinetic framework adopted 
in the quasi-linear transport model presented in this work, is believed to be a necessary 
feature for advances in first principle tokamak transport modeling. 
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Chapter 3 

The quasi-linear approximation 



In this chapter the new quasi-hnear gyro-kinetic transport model QuaLiKiz is presented. 
More in particular the attention will be focused on the quasi-linear approximation used for 
evaluating the plasma response for the transported quantities. In fact, the structure of the 
quasi-linear turbulent flux results composed by two main parts, according to the scheme: 

QL flux oc QL response (X) saturated potential (3-1) 

The first term represents the quasi-linear plasma response while the second one is the inten- 
sity of the saturated fluctuating potential. 

In the first section, the equations allowing to derive a quasi-linear response for the trans- 
ported quantities, are introduced. Particular attention will be devoted to the ordering of the 
characteristic times and to the presence of kinetic resonances. On one hand, the ordering of 
the typical times governing the plasma dynamics, defines the limits where the quasi-linear 
approximation is expected to be valid. On the other hand, the contribution from both res- 
onant and non-resonant terms appears to be crucial in the quasi-linear expression for the 
tokamak turbulent flux. 

Once the quasi-linear formulation has been established, the hypothesis of a linear response 
for the transported quantities (particles and energy) to the fluctuating potential has to be 
validated against comprehensive nonlinear simulations. This challenging aspect has been here 
treated according to different levels of analysis, namely the analysis on: (1) the characteristic 
quasi-linear and nonlinear evolution times, (2) the phase relation between the fluctuating 
potential and the transported quantity, (3) the overall transport weights (independent of the 
fluctuating potential intensity) in both nonlinear and quasi-linear regimes. 

3.1 The quasi-linear energy and particle fluxes 
3.1.1 Quasi-linear ordering and hypotheses 

After more than 40 years from the first pioneering papers [98l [29] , quasi-linear theory ( QLT) 
remains still nowadays an open subject of research that can provide a very powerful instru- 
ment for plasma physics understanding. Very large amount of literature often accompanied 
by controversies on the formulation and on the limits of QLT has been produced. Significant 
reviews can be found for example in [67l [64l [33l [7j . Even if most part of the theoretical efforts 
in QLT has been applied to ID plasma turbulence, several QL transport models have been 
proposed for the tokamak relevant 3D drift wave turbulence, providing feasible and com- 
monly used predictive tools for the evolution of the thermodynamic quantities in tokamak 
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plasmas. Among them we recall here GLF23 [M], TGLF [22], IFS-PPPL [SI], MMM95 [lU], 
Weiland model [104j . Such efforts are nowadays more and more supported because, despite 
the apparently crude approximations adopted, QLT has revealed for a relevant number of 
cases a good agreement with both experimental results [611 llOSj and nonlinear gyrokinetic 
simulations [55 ] lOTl [92] [70] . 

Inferring from the linear structure of the turbulence significant information on the nonlinear 
regime represents a great challenge; hence extremely accurate care has to be used in under- 
standing the hypotheses of the underlying physics. 



The general framework of QLT is a mean field theory for the evolution of a time averaged 
distribution function /o = {f)^ of the plasma particles population, governed by a Vlasov 
equation (2.28). A first time ordering concerning the time scale used for the distribution 
function averaging is the following: 

• r > 1/7, where 7 is the dominant unstable mode growth rate. In other words, the 
characteristic time scale t is larger than the typical time of the microscopic fluctuations. 

• T < To, where Tq is the equilibrium evolution time. This means that the characteristic 
time scale r is smaller than the macroscopic equilibrium evolution time. 

A Fourier perturbative approach on both the Hamiltonian H ~ Hq + 6h and the distribution 
function f = fo + 5f (introduced by Eq. ( 2.40 1-( 2.41 1), is apphed assuming: 

&f 6h 
— « — < 1 
/o ^0 

Averaging the Vlasov equation over the time r leads to: 



dfo 

dt 



(3.2) 



(3.3) 



while in the 5f equation, when neglecting the terms that are quadratic in the fluctuations 
amplitude, the following expression is obtained: 



dt 



(3.4) 



Using the angle-action variables (a, J), the nonlinear terms of (3.3) can be written as: 

dt /-^^ d^c 





^ dt 

T 



(2^) 



(dc^shdjSf - djShdc^Sf) 



2tt 



d^a 
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SfdaSh 



(3.5) 



QLT means that the fluctuating part of the distribution function appearing in Eq. (3.5) 
is approximated by a linear coherent response, i.e. using the important Eq. (2.44 1. Hence, a 
fundamental quasi-linear diffusion equation can be derived: 



dfo 
dt 



= V 



J • 



n Im 



n ■ djfo 
n • r2oj + ^0^ 



\SK 



(3.6) 



In the Eq. (3.6), the single-particle propagator Gn,ui = 1/ (w — n • Ooj), dependent on both 
the wave number and the frequency, is of central importance. The appearing frequency lu 
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is in principle real, nevertheless u) = uj + iO"'" is introduced in the particle propagator. The 
strictly positive additional term +iO^ represents the imaginary quantity following the Landau 
prescription on the resonance, needed for granting causality in the process. With this aim, 
since the temporal structure of the fluctuations is in the form e*"* , it is sufficient adding to the 
real frequency lo a strictly positive but infinitesimal growth rate +^0"*" , such that fluctuations 
cancel out for t — > — oo. On the other side, in the limit of 0+ 0, the term Im {G-n^^) gives 
origin to a singular Dirac function; thus meaning that in presence of a discrete spectrum 



in n and to the quasi-linear diffusion coefficient of Eq. (3.6) is not properly defined. This 
aspect is often arbitrarily solved passing from a discrete summation over n to an integral 
over a continuous space, thus recovering a proper mathematical expression. Nevertheless, 
since diffusion results from an intrinsically irreversible process, the appearance of a diffusion 
coefficient from the QL equations (which are a priori symmetric by time-reversal) is not a 
trivial point and deserves additional care. 

Substituting the real frequency w with w in the particle propagator does not simply fulfill 
causality through a vanishing imaginary contribution. Conversely, the diffusion coefficient of 



Eq. (3.6 1 has to be necessarily linked to intrinsic nonlinear effects leading to irreversibility 



through stochastic mixing of the particles orbits in the phase space. The hypothesis 

w = cj + iO+^a; + w v>Q (3.7) 

where v is definitely non negligible and positive, is actually the key point for passing from a 
conventional resonance localized QLT to a renormalized QLT. 

The stochastic character of this renormalized QLT, which results into a turbulent particles 
orbit diffusion, can be made clearer in terms of a Chirikov criterion |24) for a chaotic state; 
that is simply expressed by: 

a = I > 1 (3.8) 
where 5 and A are respectively the width of the resonant islands and the distance between 



them in the phase space. In terms of the Chirikov criterion (3.8), the non-resonant contribu- 
tions to QL diffusivity arise from the distortions of Kolmogorov-Arnold-Moser (KAM) torus 
in the phase space: a regime where the particles trajectories are not anymore integrable is 



gained when condition (3.8) is fulfilled and the pathological resonance localized character of 
diffusivity is broken in favor of a large scale diffusive behavior. Since 5 oc the resonance 

overlapping can take place either when resonances are getting closer, or when the amplitude 
of the fluctuating potential is raising. One would however expect that when moving towards 
strong turbulence with increasing \8h\, a limit in the amplitude of \5h\ exists beyond which 
QLT is not valid anymore, because of the previous assumption of 5h/ Hq <^ 1 and of the 
linearization of the particles trajectories. Nevertheless, this seems still to remain an open 
point of QLT, in particular for 3D drift wave turbulence, which could be investigated by 
future full-/ gyrokinetic simulations. A proof of the validity of QL equations in the strongly 
nonlinear regime has been indeed given within the framework of Hamiltonian dynamics in 
[34j for ID Langmuir turbulence. 

Within a renormalized QLT then, accounting for iv inside the particle propagator G^ ,^ 
effectively corresponds to a stochastic renormalization to nonlinear effects. Historically, this 
has been at the origin of the so called resonance broadening theory RBT, firstly initiated by 
Dupree in [30] and followed by several other works, leading also to more elaborate theories 
like the direct interaction approximation DIA [80l [SI] I107j . The reason of the name reso- 



nance broadening can be clearly understood since the term Im(G'„.Lj) of Eq. (3.6 1 produces 



a broader resonance function with non negligible in contrast with the singular resonance 
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localized expression found for 10+ — > 0: 



Im(G„,^)-- "^^TrJl^-n-J^oj) (3.9) 

[Lu-n-iloj) 

where the limit passage from a Lorentzian broadening around the resonance to a local 5- 
function is shown. In other words, the RBT accounts for both the resonant and non-resonant 
contributions to QL diffusivity, while the limit iO"*" — includes only the resonant terms. 

A fundamental question regards the value of such nonlinear resonance broadening v. In 
the current framework of renormalized QLT, a rather crude nonlinear dispersion relation is 
considered by: 

L" + i(7i"-^n) (3.10) 



= CJ, 



Eq. (3.10) describes a simple dynamics, where, for each mode, the linear growth rate 7^*" 
competes with the nonlinear damping i^n. Within this picture, saturation is then linked 
to the resonance broadening mechanism, when the nonlinear dissipation balances the linear 
instability drive, giving 

^QL,n = 7n" (3.11) 



Relation (3.11) actually represents the renormalized QL practical recipe for accounting the 



resonance broadening in Eq. (3.9). This kind of assumption is in fact widely used inside 
several quasi-linear transport models such as GLF23, TGLF, Weiland model, IFS-PPPL, 
MMM95. 

The stochastic nature of the present QLT can be regarded in terms of a random walk 
phenomenology. Within this framework, one can estimate the perpendicular diffusion (cross 
field transport in the case of tokamak plasmas) as: 

D^^\i^\5vi_\^)T^j: (3.12) 

where 5v±^ are the velocity fluctuations in the perpendicular direction (the E x B drifts in 
tokamak plasmas); r^,p is a characteristic wave-particle interaction time, or more intuitively 
the effective lifetime of the field pattern. In order to test the validity of such QL model, one 
should compare this time scale r^p with a characteristic nonlinear time t^^. QLT applies 
in fact when the particles motion perturbations remain small during r^p, while one expects 
that for times t > tnl the linearization of the trajectories and the QL diffusivities will not 
be valid anymore. This means that a safe condition of validity for this QL random walk 
representation results from the following ordering of these time scales: 

Twp < tnl (3.13) 

The latter criterion reflects the need of avoiding that particles undergo to a trapping condition 
in the resulting electric field pattern, which is incompatible within the present QL formulation 
It is however important to notice that also different approaches accounting for turbulent 
trapping have been developed for example with the so called clump theory [3^ |5J [M] . 

The characteristic nonlinear time can be estimated starting from a general definition of 
the diffusion coefficient in the velocity space [Ml [53] : 

D = lim \ ^ - / C (r) dr (3.14) 



^The use of the word trapping, originally introduced by Dupree, could however be misleading, since 
coherent phase space islands are destroyed when considering a spectrum of waves. 
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where Vq is the velocity initial condition and the average (. . .) is to be taken over an ensemble 
of particles. The appearing correlation function is defined as: 



C{t) 



(5E [x {to + r) , to + r] [x (to) , to]) 



(3.15) 



Relation (3.151 actually defines a Lagrangian, taken along the orbits, correlation function 



for the fluctuating electric field i5E. These correlations will decay on a nonlinear time scale 
Ti , measuring the decorrelation time of the turbulent structures experienced in a framework 
moving with the particles velocity. Indeed, one can demonstrate [63j that the time integration 
according to (3.141 lead to a velocity diffusion coefficient similar to the one embedded in Eq. 



(3.6 1, with a resonance broadened function for the single particle propagator (in general not 
necessarily with a Lorentzian shape). The principal effect of this kind of approach is then the 
introduction of a characteristic resonance broadening width = t'l 1 ii^ the single particle 
propagator Gn.u- 

The former discussion about stochastic QLT can then be reviewed in the light of the latter 
implications. In particular, a useful criterion for the validity of a simple QL random walk 



model can be gained if testing the QL ordering (3.131, where the nonlinear time coincides 
with the just defined Lagrangian one: t^l = tl- In the case of tokamak microturbulence, 
this kind of approach has been used for example in |70j . by mean of massive nonlinear 
gyrokinetic particle in cell (PIC) simulations of ETG turbulence. The same argument will 
be also deepened later in this work. 

Another kind of approach can rely on the comparison between the resonance broadening 
width adopted in the renormalized QL model according to Eq. (3.11) with the nonlinear 
fNL- That would imply that a nonlinear Lagrangian time is expected to be close to the 
inverse of a linear growth rate: 



tnl- 



1 

In 



(3.16) 



However this point should deserve great additional care for a further improvement of QL 
transport models, and represents an open issue also discussed in the following. 

Heuristic derivation of a saturation rule 



Until here the main hypothesis of the present QL formulation has been the stochasticity 
of the particles trajectories in the phase space, thus allowing a random walk model. The 
second main critical issue for gaining QL predictions for the turbulent fluxes is the saturation 
level of the fluctuating potential. This argument will be the subject examined in Chapter 
4. Nevertheless, in the present framework, it is worth noting that a heuristic derivation of 
a commonly adopted saturation model, namely the mixing length rule, can be coherently 
derived from the latter considerations. 

Following from the former discussion on the QL hypothesis of resonance broadening driven 
saturation Eq. (3.10)-(3.11 ), one expect that also the saturation level is linked with the QL 
resonance broadening width vql.h- At this regard a significant conclusion can be gained if 
supposing, without a rigorous justification, the equality between i'QL,n and the inverse of 
nonlinear decay time of the Lagrangian correlation function (3.15). An estimation of this 



time scale can be obtained making use of a quite strong hypothesis of random Gaussian 
statistics for the saturated fluctuating electric field [64] : 



1 

TNL,k 



(kl) 



(3.17) 
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where the nonhnear decay time scales like the inverse of a Dupree-Kolmogorov-like time 
[231 [M] • It has to be stressed that the assumption of Gaussian statistics leading to relation 



(3.17) constitutes a strong additional hypothesis, consistent with a Markovian limit of a 
diffusion equation. A priori, this kind of choice has not any rigorous justification, since for 
developed plasma turbulence the fluctuating potential will be self-consistent with the particle 
distribution function through the Poisson- Ampere equation. Finally it is just sufficient of 
using the already discussed equality VQL,n ~ ^/tnl.ii to obtain: 



(3.18) 



Eq. (3.181 is usually referred as a mixing length saturation rule: within this renormalized 



QLT, it arises from the resonance broadening. The rule (3.18) is a very important practical 



relation adopted in several QL transport model; despite it is not rigorous, there are several 
indications showing that this kind of hypothesis can at least qualitatively reproduce nonlinear 
fluxes. Recently, a remarkable confirmation of this kind of behavior has been identified also 
by mean of comprehensive nonlinear gyrokinetic simulations of tokamak TEM turbulence 

m. 



From expression (3.181 it also clearly appears that the hypotheses adopted in deriving the 
estimate (3.17) are crucial. Non-Gaussian deviations of the saturated field statistics could in 



fact result in a perpendicular diffusion coefficient D± which deviates from the linear scaling 
with 7^™, as prescribed by the usual mixing length estimate (3.18). Finally one should 



notice that the arbitrary dependence of the QLT on the wave-vector n is here appearing in 



expressions (3.17)-(3.18), but this argument will be deepened in the following. 



3.1.2 A new quasi-linear gyrokinetic model, QuaLiKiz 

The quasi-linear framework just described has been applied to the calculation of the turbulent 
fluxes of energy and particle in tokamak plasmas. The gyrokinetic linear solver Kinzero ]12], 



introduced in the previous paragraph 2.2.2 has been upgraded to the quasi-linear transport 
model named QuaLiKiz [22] [13] . 

The gyrokinetic quasi-linear formulation of the particle F^ and energy Qs flux for each 
plasma species s, follows from the velocity moments integration of the diffusion coefficient 



embedded in Eq. (3.6). The derived expressions are: 



Fs = Re (^5ns ^^^^^ 



3\ RVrTs 



Im 



(3.19) 
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The sum is over the toroidal wave-numbers n, while the frequencies are: 

k\\v\\ w ^ —wVt kg = — («J-21) 

" " q R r 

The details on the notation and the integration over the passing and trapped domains are 
identical to the ones introduced for the description of Kinezero in the paragraph |2.2.2[ Thanks 
to the relation kg ~ nq/r, there is a correspondence between the toroidal wave-numbers 
labeled by n and k, which will be from here equivalently used in the text. 

A first relevant consideration stems from the unphysical breakdown of the particle flux 
ambipolarity when using a strict resonance localized quasi-linear diffusion. Considering for 
simplicity a plasma of ions and electrons only, the quasi-neutrality condition implies Srie^k — 
Srii^k for a given wave-number k: thus automatically turns into the ambipolar fluxes Fg = 



Ti^k- Making use of the useful resonance limit (3.9), for 0+ -> 0, the quasi-linear expression 



(3.19) implies for a given k: 



r (3.22) 



^nf2s {£, A) 

nujds J ' ' / £.x 

Since the ion and electron resonances do not normally overlap, it appears that the particle 
flow is not automatically ambipolar, i.e. Tf^ ^ 7^ Tj.fe- This breakdown following from Eq. 



(3.22), makes that quasi-linear resonance localized formulation is usually rejected in the 



tokamak transport modeling in favor of a renormalized QLT, introduced in paragraph |3TTT 



In principle, two kinds of broadening exist in the quasi-linear fluxes expressed by Eq. 



(3.19)-(3.20), whose physical origin is distinct and has to be specified. The first one actually 



coincides with the just mentioned RBT (see Eq. (3.9)): as already explained, it corresponds 
to a addition of a non-negligible term -t-iO"'' — > in the particle propagator Gn,u:- The 
second one is instead related to an intrinsic w-spectral shape of the fluctuating potential 
|^<?!'n,cj|'- Since the properties of |(5(^„,cj|' derive from the nonlinear saturation, it is reasonable 
to assume that this spectral quantity can deviate from a singular Dirac function localized 
at the frequency of the linearly unstable mode; conversely, it can exhibit a finite broadening 
around a given frequency. Here we refer to this second broadening mechanism as frequency 
broadening, in order to discriminate it from the previous resonance broadening. General 
models on the |5</)„,aj|' frequency spectrum are commonly considered by statistical plasma 
turbulence theories [64l I102| : a relevant example is the so called pole approximation of the 
Lorentz line broadening model used by the DIA theory [56l I102j . Letting wofc be the linear 
/c-mode frequency, this model takes the following form: 

\5^kJ = \5^u\' ^ (3.23) 

where wq/c + presents a nonlinear shift with respect to the corresponding linear frequency, 
while Qffe is the frequency broadening or the nonlinear decorrelation rate. 

A rather general approach that can be adopted in tokamak transport quasi-linear models, 
is describing the \5(j)k,uj\^ spectrum through a sum of w-broadenings around the different linear 
eigenmodes labeled by j: 

\5<t>kJ = \5<t>k?Y.^k,, (^) Sk, {uj) = (3.24) 

(w-wofej) +ai. 
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It is worth noting that with respect to the DIA model (Eq. (3.23)), the nonhnear shift 5uik 
has been skipped. CoupHng in the quasi-Unear particle flux (3.19 1 both the frequency (Eq. 
(3.24)) and the resonance broadenings (Eq. (3.9)) would end up in the following expression: 



r,, (X 



E 

n 

E 



7 

duj 

— 1 

TT 



e^^Im 5fej ( 



Lo - nfls {£, A) + iO+ 
Lu — nuj 



uj — nQs (^! A) + ii'kj 



e.x 



1 



(3.25) 

(f — I) ; the integration over 



where for simplicity we have used nuj* = ^ ^ 

duj can in principle be analytically performed by residues. 

At this regard, most QL tokamak transport models assume more or less implicitly that 
for each wave number k exists a well defined frequency lo such that to — > ujq}~ . In other words 
this choice corresponds to: 

(w) ^ (5 (w - wofcj ) (3.26) 

On the contrary, QuaLiKiz explicitly assumes the Lorentzian frequency broadening expressed 
by Eq. (3.24), while keeping the localized resonance limit zO+ — > 0. The combination of 
these choices results in a formulation completely equivalent to the more common renormal- 
ized quasi-linear theory, i.e. properly accounting for both the resonant and non-resonant 
contributions to the quasi-linear flux. Indeed, it is possible to write: 



E 

°^"-E 



\f£e ^Im ( Sk,j {uj) 



Jo {k±Ps 



"fe, 



nils + *0 



|<5</>„.jr (3.27) 
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Jo' [k^ps) 1 50, 



(3.28) 

The latter demonstration proves finally that the frequency broadening model proposed in 
QuaLiKiz (Eq. (3.27)), coincides with the renormalized QLT (Eq. (3.28)) used in several 
other tokamak quasi-linear models like GLF23, TGLF, Weiland model, etc., under the fol- 
lowing conditions: 

1. Sk.j (i^) is chosen with a Lorentzian shape. Any different functional dependence for 
Sk,j {^) would have in fact lead to a discrepancy between the two approaches. 

2. oikj — i^kj- The broadening introduced in the renormalized QLT through the prescrip- 
tion Lj — !■ LjQkj + ivkj is identical to the Lorentzian frequency spectrum width oikj ■ 

At the actual time, neither the formulation proposed in QuaLiKiz, nor the more common 
renormalized QL models, simultaneously accounts for both these broadenings, as suggested 
by Eq. (3.25). In principle both these widths, namely the resonance and the frequency 
broadenings, could be object of study through dedicated nonlinear gyrokinetic simulations. 
In practice, separating these two effects is not an easy task, because the fluctuations of the 
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potential and of the particle distribution function are back-reacting one on the other through 
the Poisson- Ampere equation. 

A crucial point for all the actual quasi-linear transport models is the choice for the broad- 
ening Ukj- The practical solution involves the balancing of the linear growth rate, relying on 
considerations expressed by Eq. (3.111, then giving: 



Q-kj — Ikj 



(3.29) 



Eq. (3.29) is presently adopted by QuaLiKiz as well by all the other tokamak quasi-linear 



transport models. Nevertheless, this hypothesis represents an assumption that is not rigor- 



ously justified. A remarkable example where the rule (3.29) breaks down can be easily found 
in the case of a linearly stable mode at a given wave-number k, i.e. 7^ — 0. Whereas in 



this case quasi- linear models find = following from (3.29), nonlinear simulations |100j 



show instead a finite value of the turbulent flux ^ 0. One way to account for these kind 



of non-local couplings in the k space would be to replace the rule (3.29 1 by the non-local one 
Oikj = Ikj +k^D. In that way, even if 7^, = 0, a finite contribution to the quasi-linear flux will 
be retained. The additional term +k^D, where _D is a general total diffusivity independent of 
k, can still be heuristically justifled through statistical plasma theories [64]; more speciflcally 
the nonlinear time derivative becomes d/dt — > d/dt — k^D, resulting in the modifled linear 
growth rate 7^, = 7;^ -|- k^D. 

On the other hand, also turbulence measurements can provide useful insights about this is- 
sue, and the argument will be deepened in Chapter [4j 



3.2 Validating the quasi-linear response 

The purpose of this section is the validation of the assumption of the linear response of the 
transported quantities to the fluctuating potential. The only way to quantify the goodness 
of the quasi-linear response for the turbulent energy and particle transport is an accurate 
comparison with comprehensive nonlinear simulations. 

Despite the large number of existing tokamak quasi-linear transport models, most of them 
are commonly applied, more or less successfully, for predictive simulations of plasma the dis- 
charges (i.e. calculating the time evolution of the Te, Ti and profiles), without a clear 
picture of the capabilities and the limits of this kind of modeling. 

Only very recently, for example in the case of the TGLF model, an exhaustive comparison of 
the total turbulent fluxes has been performed between the quasi-linear (TGLF) predictions 
and the full nonlinear gyrokinetic expectations |59) . Even if this work represents a great 
advance, the relevance of this validation mostly relies on the specific choices adopted in the 
quasi- linear model that contribute to the final estimation of the total flux. In fact, when a 
theory based model fails to match the nonlinear gyrokinetic simulated transport (or the ex- 
perimental transport for that matter) , it is generally not possible to know if there is a failure 
in the accuracy of the underlying physics (e.g. gyrofluid versus gyrokinetic description), of 
the nonlinear saturation rule, or of the quasi-linear response. 

In this thesis work, the validity of a tokamak relevant quasi-linear approximation is sys- 
tematically studied apart from any hypothesis on the nonlinear saturation mechanism. By 
using both linear and nonlinear gyrokinetic simulations, the aim is then to isolate and quan- 
tify the success or failure of the quasi-linear approximation itself. This analysis is believed to 
be a crucial point to get insight about the effective predictive capabilities of the model. The 
validation of the quasi-linear response is structured according to the following arguments: 
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• In order to verify if the quasi-linear expected time ordering is respected, the charac- 
teristic nonlinear and linear times are defined and compared in nonlinear gyrokinetic 
simulations of tokamak microturbulence . 

• The cross-phase between the fluctuating transported quantity (energy and particle) and 
potential are compared in the nonlinear and in the linear regime. 

• An overall transport weight is defined and a relevant ratio between the quasi-linear and 
the nonlinear response is quantified. 

The reference case for the nonlinear simulations 



The most part of the nonlinear gyrokinetic simulations presented in this work are per- 
formed using the code GYRO 19j. Further details about the GYRO code and the nonlinear 
gyrokinetic simulations are reported in Appendix [B] A standard reference case of study is 
defined in Table 3.1 and it will be later referred as GA-standard case] here below a summary 
of the physical parameters which are adopted; The numerical grid and the model assump- 
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Table 3.1: Plasma parameters defining the GA-standard case. Unless otherwise specified, 
this set refers to electrostatic turbulence in circular s — a magnetic geometry. Moreover, 
ctMHD = 0.0 and Z^ff — 1.0. 

tions for the GYRO nonlinear simulations are here briefly summarized (otherwise specified, 
these parameters apply to all the GYRO simulations presented in this chapter): 

• Local (flux-tube) simulations with periodic radial boundary conditions 

• Drift-kinetic electrons (electron FLR are neglected) with real mass ratio y^rrii/ nie — 60, 
electrostatic, coUisionless 

• Box size in the perpendicular directions [L^/ps, Ly/ ps] — [126,126], radial resolution 
Ax/ps = 0.75 

• 16 Complex toroidal harmonics, covering 0.0 < kyPs < 0.75, 12 grid points in the 
parallel direction, 15 in the gyroaverage and 5 in the radial derivative 

• 128-point velocity space discretization per spatial cell, 8 pitch angles, 8 energies and 2 
signs of velocity 

• Statistical averages refer to typical time intervals 100 < t < 1000, where time is ex- 
pressed in a/cg units 

In the GYRO code, upwind dissipative advection schemes are used in order to provide the 
dissipation and time irreversibilty required for the achievement of statistically steady states of 
turbulence; this numerical dissipation occurs only in the real space and arises from the upwind 
operators. Hence GYRO does not use any velocity-space dissipation other than possibly 
the collision operator. It has been shown |20| that adding upwind dissipation to radial 
advection terms is required to smooth over sub-grid-scale numerical disturbances associated 
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with electron Landau layer physics. The radial upwind differencing in the drift-advection 
terms is then necessary when solving the electron equations on a p^-scale grid. In cases 
where electron Landau layer effects do not need to be resolved, upwind dissipation eliminates 
unwanted numerical effects which would occur in a non-dissipative schemes. 

3.2.1 Characteristic turbulence times 

Historically, the quasi-linear theory has been elaborated for test particles [M]. As already 
pointed in paragraph 3.1.H this point can be easily understood, since the quasi-linear theory 



does not provide a self-consistent treatment; this is due to the fact that there is no back- 
reaction of the perturbed quantities on the fluctuating potential. 

A powerful tool derives from fluid turbulence theories, the so called Kubo number JC [|6&, a 
property that can be defined for an advecting velocity field. The Kubo number is the ratio 
between (1) the wave-particle interaction time T^p characterizing the lifetime of the pattern 
the particle senses, and (2) a flight time T/(or eddy turnover time), characterizing the time 
a particle would spend around the field structure, t/ can also be interpreted as the time for 
the macroscopic flow to advect a perturbation across the system. Then, 

_ wave - particle time 

/t = — : = (o.oUj 

eddy turnover time tj 

The main difficulty comes from the fact that usually the advecting field is a function of both 
the 3D-position and time, so that r^p should be a Lagrangian time, which is typically hard 
to compute. According to the statistical plasma turbulence theories, a condition of validity 
for this quasi-linear framework is that the particles should not be trapped in the field (see 



paragraph 3.1.1 Eq. (3.8) and Eq. (3.13)). In terms of the Kubo number, this condition 
turns in: 

K. < 1 (3.31) 



Note that quasi-linear ordering set by Eq. (3.31) through the Kubo number is completely 



analogous to the condition expressed by Eq. (3.13). It is worth noting that there are examples 



where quasi- linear theory has been shown to work up to /C « 1 |81j . Ideally one would 
compute this Kubo number from nonlinear gyrokinetic simulations. 

Since the tokamak nonlinear gyrokinetic simulations consistently compute the coupled 
Vlasov-Maxwell equations, the definitions of these characteristic times have to be reviewed. 
Tf can be calculated as the ratio between an auto-correlation length L^, and the velocity at 
which the fluctuating quantities are radially transported Svr'. 

Referring to the tokamak geometry, Lc is evaluated as a radial correlation decay length. The 
definition of the latter one requires the introduction of the 2D correlation function: 



' (r, (fi, t) 6(f>* (r + Ar,(f + Atp, t))^ 

|2\ 



a.^ (Ar, A<^) = """'"'7" " ' '7k ' '^'^ (3-33) 



I r.ip.t 

The radial correlation function Cr (Ar) is calculated by taking the maximal value along the 
ridge of Cr.ip.t (Ar, A(^). The function Cr (Ar) typically presents an exponential decay over 
the radial separation Ar, as shown in Fig. |3.1[ 
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Figure 3.1: Example of a radial correlation function C,. (Ar) obtained from a nonlinear 
GYRO simulation using kinetic electrons and collisions. The dotted line is the exponential 
fit. 
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Finally, the radial correlation length is defined as = i^I^iiiii. 

Tttip, the effective lifetime of the field pattern, is estimated according to the arguments leading 



to Eq. (3.12): together with the expression of the radial velocity fluctuations at a given 



toroidal wave-number fc, can be written as: 
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(3.34) 



where D is the particle diffusivity, D = — Fe ^/VrJi. This approach has also been applied in 
Ref. [TD], in the case of nonlinear simulations of ETG tokamak turbulence. 

In the present work, this Kubo-like number defined by Eq. (3.32) is computed from non- 



linear gyrokinetic simulations of coupled ITG-TEM turbulence using GYRO. With reference 
to the GA-standard case, a wide scan on the normalized temperature gradient R/Lt has 
been explored, varying simultaneously R/Lxi = R/Lxe- The Kubo-like number /C computed 
from each GYRO simulation is shown in Fig. |3.2[ For these parameters, /C exhibit val- 
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Figure 3.2: Kubo-like number /C (Eq. (3.32)) computed from GYRO nonlinear simulations 
versus a R/Lt scan on the GA-standard case; vertical bars in the graph refer to the statistical 
evaluation of the intrinsic turbulence intermittency predicted by the nonlinear simulations. 



ues well below 1.0 even for highest gradients in the scan. These estimations provide then 
a useful information on the expected quasi-linearity character of tokamak relevant micro- 
turbulence, involving both ion and electron unstable modes. Interestingly, the quasi-linear 
criterion K, < 1 (Eq. (3.31)) appears here to be largely respected. The numerical solution 
of the full nonlinear dynamics within the gyrokinetic framework has then revealed that the 
particles do not undergo to a trapping condition in the fluctuating electric field pattern. 
Hence, in the parametric region here explored, these results justify further attempts towards 
the quasi-linear modeling of the tokamak turbulent fluxes. 
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3.2.2 Nonlinear versus linear phase relations 

Apart from the criterion based on the ordering of the turbulence characteristic times just 
described, another relevant test for verifying the hypothesis of the quasi-linear response relies 
on the phase relations between the fluctuating fields. Referring for example to the particle 
flux, the general nonlinear expression for a given wave-number k can be written as: 

rfc=Re((5nfc — 'b I ^ V — B — smA<I>^ (3.35) 

where denotes the cross-phase between the fluctuating particle and potential fields. 



Expressions analogous to Eq. (3.351 can be written for both the energy fluxes, defining the 



ion and electron energy cross-phases, respectively A$f^' ^'^ and A<i>f^° The representa- 



tion of Eq. (3.35) is particularly physically meaningful, since it highlights that the turbulent 
transport can only arise in presence of phase shift between the fluctuating potential and the 
transported quantities. Phase shifts (or cross-phases equivalently) close to or tt are then 
indicative of a small fluxEl 

The quasi-linear theory can not give any insight on the product of the fluctuations abso- 
lute values l^rifcl while it provides well defined phase relations between the fluctuating 



fields (see Eqs. (3.19 l-(3.20|j^| The comparison of the cross-phases in the linear and the 



nonlinear regimes allows to verify the validity of the quasi-linear response hypothesis. Nev- 
ertheless, as it will be deepened in the following, the linear cross-phases represent a relevant, 
but still not the whole information carried by the quasi-linear approximation. 

Practically, this aspect is investigated through both nonlinear and linear gyrokinetic sim- 
ulations using the GYRO code. In the nonlinear regime, the cross-phase is not a constant 
quantity, but it presents oscillations reflecting the statistical intermittency of the turbulence. 
For this reason, the nonlinear cross-phases A$^' will be presented as probability density 
functions (PDF) of the phase shift between the fluctuating quantities computed by the sim- 
ulation, i.e. between [(Srifc (r, £) , 5Ei_}^ (r, i) , SE^^k {r, 0, t)] and (5(/>^ (r, 9, t) ; where ^ = 0, 
r spans over the whole radial domain and a nonlinear saturation t-interval is considered. 
Conversely, only the linear most unstable mode is used for the estimation of the quasi-linear 
cross-phases, which will be referred as A$|'-^. 



Fig. [3^ presents the results from the GYRO simulations on the GA-ITG-TEM stan- 
dard case, which nonlinearly finds dominant (at low fcgPs < 0.4) saturated modes in the ion 
diamagnetic direction. Previous works in the plasma literature have already analyzed this 
issue for pure electron TEM turbulence [371 [SS]. It appears that the linear phase shifts, i.e. 



the white lines in Fig. 3.3 remain close to the maximal PDF values of the nonlinear cross- 
phases. A very good agreement is therefore simultaneously observed for coupled ITG-TEM 
turbulence for all the particle and energy transport channels. It has to be however noticed 
that the quasi-linear phase shifts significantly breakdown at high k. Nevertheless, this failure 
is expected to have a small impact on the goodness of the quasi-linear approximation. As 
quantitatively detailed in the following in fact, most part of the transport is driven at scales 
corresponding to k±ps ~ 0.2, making that the cross-phases for k±ps > 0.5 contribute less to 
the total turbulent flux. 

It is important to check that the agreement on the cross-phase relations between the lin- 
ear and the nonlinear regimes is respected across a variation of the plasma parameters. Hence, 



two modified GA-ITG-TEM standard cases are analyzed. The first one, reported in Fig. 3.4 



^It immediately appears that within the approximation of adiabatic electrons, where the density response 
is simply taken as 5n/n = e5(f>/Te without any phase shift, no particle flux can be retained. 

^The linear dispersion relation used to derive the quasi-linear fluxes naturally implies the definition of a 
complex quantity, with a well defined phase relation between S<l> and [Sn, SEi, SEe] 
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Figure 3.3: PDF of the nonhnear cross-phases (color contour plot) and the linear cross-phase 
of the most unstable mode (white line): a) dn — Scj), b) 6Ei — S(f> and c) SE^ — S4> from a local 
GYRO simulation on the GA-ITG-TEM standard case. 
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Figure 3.4: PDF of the nonlinear cross-phases (color contour plot) and the linear cross-phase 
of the most unstable mode (white line): a) Sn — S(p, b) 6Ei — 5(f> and c) 6 Eg — S4> from a local 
GYRO simulation on the modified GA-ITG-TEM case with R/Lrt.e = 9.0 ^ 18.0. 



50 



The quasi-linear approximation 



5n-6{t) 



5E.-5^ 




-10 12 
angle [rad] 



3 _3 _2 -1 
angle [rad] 



1 -3 




-2 -1 
angle [rad] 



Figure 3.5: PDF of the nonhnear cross-phases and the hnear cross-phase of the most unstable 
mode (white hne): a) Sn — Scj), b) SEi — Scj) and c) SE^ — 5<t> from a local GYRO simulation 
on the modified GA-ITG-TEM case with R/Lti = 9.0 ^ 6.0. 



is intended to test the quasi-linear response when enhancing the turbulent drives by increas- 
ing both the ion and electron temperature gradients, i.e. R/Lxi = R/Lxe = 9.0 ^ 18.0. 
As shown in Fig. |3.4[ the linear phase shifts relative to the electron and ion energy trans- 
port appear still reasonably close to the nonlinear values, while a more pronounced failure 
is found for the particle transport. With these parameters, the GYRO nonlinear simulation 
predict a strong inward flow, mostly due the dominant low-fc ITG turbulence. On the other 
hand, since the values of the linear dn — Scj) cross-phases are more positive than the nonlinear 
phase shifts, the quasi-linear estimate for the particle flux is possibly expected to disagree 
with the nonlinear flux for a non negligible over-prediction. The quasi-linear particle flux, 
especially for strong turbulence cases and inward flows, deserves then a careful treatment: a 
more exhaustive and quantitative approach will be introduced in the following paragraph. 
A second case explore turbulence conditions where the ion ITG turbulence drive is lowered. 
This is realized changing only the ion temperature gradient with respect to the GA-ITG- 
TEM standard case, i.e. R/Lti = 9.0 6.0, where ITG and TEM unstable branches are 
in competition with similar linear growth rates: results are shown in Fig. 3.5 Here the 



breakdown of the linear cross-phases with respect to the nonlinear phase shifts is even more 
evident. The most relevant failure is again on the particle transport, where the linear result 
significantly departs from the nonlinear PDF. The jumps observed in the linear Sn — dcf) cross- 
phase are due to the variations of the linear most unstable mode between ion and electron 
diamagnetic directions across the kg wave-numbers. Moreover, also the phase shifts relative 
to the electron transport exhibit a relevant disagreement, with the values of the linear cross- 
phases well below to the nonlinear ones. Also in this case then, these results could suggest 
that the quasi-linear estimates on the turbulent fluxes deviate from the nonlinear predictions. 
Nevertheless, it is important to remind that the linear cross-phases shown in Figs. |3.3[ |3.4| 



and 3.5 refer only to the linear leading mode, while in general several unstable mode in both 
the electron and the ion diamagnetic direction are active. These sub-dominant modes are 
expected to drive a non- negligible contribution to the total turbulent fluxes, especially in 
conditions like the ones reported in Fig. |3.5[ where ITG and TEM are competing at similar 
k scales. 
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In summary, the analysis of the hnear versus nonlinear cross-phases appears as a pow- 
erful tool for validating the hypothesis of the quasi-linear response. Here we have shown 
that a good agreement on the phase shifts is obtained for conditions of coupled ITG-TEM 



turbulence, simultaneously on the particle and energy transport channels (Fig. 3.3). The 
breakdown of the linear cross-phases generally observed at high k, is not expected to have a 
drastic impact on the total turbulent fluxes which are mainly driven at k±pg < 0.5. 
Nevertheless, two main drawbacks of the quasi-linear validation approach based on the anal- 
ysis of the cross-phases emerged; these are: 



The quasi-linear cross-phases shown in Figs. 3.3[ 3.4 and 3.5 are derived from the linear 



most unstable mode, while it is reasonable to expect that also sub-dominant modes play 
a role in the estimation of the total flux. Even if in principle the cross-phases relative to 
these sub-dominant linear cigcnmodes can be numerically solved, that would not bring 
a definitive information on the final goodness of the quasi-linear response. The latter 
one in fact, will result from the weighted contributions of several unstable branches. 
The cross-phases of the linear unstable modes do not carry any insight on these relative 
weights. A more clever way to account for the role of sub-dominant modes have then 
to be elaborated. 

As suggested by the previous point, the cross-phases only partially represent the infor- 
mation provided by the quasi-linear theory (at least within the formulation presented 



in Sec. 3.1 1. In order to quantitatively address the validation of the quasi-linear re- 
sponse, another strategy has to be adopted, able to capture the whole outcome of the 
quasi-linear approximation. 

Both these points motivate dedicated approaches to the validation of quasi-linear modeling 
and they will be originally treated in the next paragraph. 

3.2.3 Quasi-linear transport weights 

The nonlinear particle as well ion and electron energy fluxes [Te,Qi,Qe] and their respec- 
tive effective diffusivities [De^XiiXe] are given by the correlation of the particle and energy 
moment fluctuations per wave-number k [Sne,k,SEi^k,SEf,^k] with the radial E x B velocity 
fluctuations Svr^k- 

= Rc ^ ( [Sn-J^„6E^j,,SE:^,] <5<;^') (3.36) 
fc 

where the superscript nl refer to actual fluctuations that can be computed for example from 
a nonlinear simulation. Conversely, the quasi-linear turbulent fluxes are obtained convolving 



a quasi-linear response with a saturated potential spectral intensity (Eq. (3.19)-(3.20)). The 
quasi-linear convolution takes the general following form: 

QL-flux cx ^ QL-weight^, ® Spectral-intensity ^. .^ (3.37) 
k,j 

The first part derives from the linear correlation of the density (or energy) perturbations 
and the radial E x B velocity perturbations. The spectral intensity captures instead the 
saturated strength of the turbulence, which depends on the nonlinear coupling of the wave 
numbers. The index j spans over a discrete number of unstable modes across the continuous 
frequency w. 
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It is important at this point to detail tlie relation that links the quasi-linear transport 



weights appearing in Eq. (3.37) and the cross-phases just previously defined: 



QL-wcightj.j = 



ke \[Snkj,SEi^k,jSEe,k,j]\ 



B 



lSn.5Ei,SE^]-S<p 



(3.38) 



where the here brackets refer to a flux-surface average. Eq. ( 3.38| ) points out the reason why 
the linear cross-phases alone do not carry the whole information on the quasi-linear response, 
as mentioned in paragraph 3.2.2 In fact, the amplitude ratios |[(5nfej-, SEi^k.j^Eg^k.jW / \S4>k,j\ 
appearing in Eq. (3.381 are effectively part of the quasi-linear transport weights, defining 



linear relative amplitude weightings for each unstable mode. 

There are three main strategies that can be adopted when practically calculating these 
transport weights, and therefore the quasi-linear turbulent fluxes, according to expression 
(3.37). These are here below summarized. 



1. Eigenvalue code approach. A discrete number of linear eigenmodes are found by a 
numerical solver; for each wave-number k there are j different solutions. QL-weight^, ^ 
is then unambiguously referred to the corresponding complex linear eigenvalue wj, j = 
i^oA; j"+i7ofc j". Defining the following linear complex functional (that immediately derives 



from Eqs. (3.19) and (3.20)) 



did 



£Pe'^Sk,j (w) 



uj — nflg (£, A) -t- ivkj 



(3.39) 



where p = 1/2 when dealing with particle transport and p = 3/2 for energy transport, 
the general expression for the quasi-linear weight is straightforward: 



QL-weightfc J = lm{J'k,j) 
In this framework, the expression of the quasi-linear cross-phase is instead: 

A$?'.=arctan-^"'(-^'='^^ 



k,j 



Re(-Ffcj) 



(3.40) 



(3.41) 



(3.39) 



From Eqs 
weights, as well of the cross-phases 



(3.40) and (3.41), it appears that the actual value of the quasi-linear 
here depend on the choice on the resonance and 
frequency broadenings, respectively Vkj and Skj retained within the linear dispersion 



relation, as discussed in paragraph |3. 1.2 The QuaLiKiz model for example, calculates 
Eq. (3.39) in the limit of Vkj 0+ and Sk,j with a Lorentzian shape. For this reason, a 
systematic verification of the quasi-linear response, apart from any further hypothesis, is 
investigated through the following two points rather than with the eigenvalue approach. 



2. mQL, Leading mode approach. In the case of initial value codes (such as the linear 
version of GYRO, sec Appendix [b| for details), only the linear leading mode can be 
numerically solved; the index j is then here fixed to j = 1. In this case, the quasi-linear 
weight can be written as: 



QL-weightj, 1 = 



Re 



X T^lin r rplin 



'k,l\ 



(3.42) 
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while the cross-phase takes the form: 



^^fi — arctan 



Im 


( 


s:^lin r rplin s: jiplin 


) 


Re ( 


'^"-fe,l'°-^i,fc,l''^^e,fc,l 


) 



Im 



— arctan 



Re 



(3.43) 



where the superscript lin refers to a linear simulation. It is worth noting that both the 
eigenmodes appearing at the numerator and at the denominator of Eq. (|3.42| are ex- 



ponentially diverging in timcj^ since they are linearly unstable: the ratio between these 
two quantities is instead properly defined. Moreover, the advantage of this approach 



is that the quasi-linear weight of Eq. (3.421 does not depend on any additional choice 
on the broadenings, contrarily to the eigenvalue approach. The time evolution of a 
normal linear mode in the initial value codes is in fact solved using the quasi-neutrality 
constraint 6n(,^k,i = <5^i,fc,ij so that the quasi-linear particle flows are automatically 



ambipolar, i.e. - e,fc,i " " i,fc,r 
Nevertheless, according to this formulation, one should be able to identify and separate 
a nonlinear spectral intensity for each unstable mode j at the wave-number fc, i.e. 



Spectral-intensity 



(3.44) 



Unfortunately, this kind of information can not be inferred from nonlinear simulations, 
since the distinguishable structure of the linear eigenmodes does not survive in the 
nonlinear saturation regime. Therefore, only one fc-resolved nonlinear saturation spec- 
trum can be computed, i.e. (|<5'/'fc'|)- In order to gain an estimation of the quasi-linear 
turbulent flux according to Eq, 

proximation 







S^kll 





(3.371, the mQL approach operates through the ap- 
It is important to detail the main consequences of 

this approximation. 

Since the quasilinear weight of the leading mode is expected to be larger than that for 
the sub-dominant modes, the full spectral intensity applied only to the leading mode, 
will tend to produce a larger quasi-linear transport than it would be obtained if cor- 
rectly distributed over both leading and sub-dominant modes. In fact, the leading mode 
in each toroidal wave-number typically refers to an ion or electron directed unstable 
branch. Moreover, within each branch, both the initial value mQL and the eigenvalue 
code approaches keep only the most outward ballooning mode, which is typically the 
most unstable. Effectively, there is a continuum of lesser ballooning and less unstable 
modes with smaller quasi-linear weight and lesser but not zero spectral weight. These 
are left uncounted, while in principle they can contribute to the quasi-linear transport. 
Practically, the accuracy of the mQL approach is tested by two-step gyrokinetic simu- 
lations: firstly, a linear run is performed to get the quasi-linear weight on the lead- 
ing mode, then a nonlinear simulation on the same parameters allows to get the 
nonlinear spectral intensity as well as the actual nonlinear transport. The quasi- 
pmQL Q^Q^ Q"^Q^ 



linear fluxes 



obtained according to Eq. 



(3.42) can then be 
The ratios 



compared to the nonlinear transport in each channel and wave-number 

Y'mQL^Y"^i^QrnQL^Qni^Q,nQLjQni] (otherwise Called overages) are physically mean- 
ingful, since they are independent of the structure of the saturation spectrum that is 
fixed by the nonlinear simulation. 



■'The linear eigenmodc temporal dependence follows Scpi^ (t) = Scfi^^ '"fc.i'+Tfc.i* 
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3. fQL, Full frequency spectrum approach. In contrast to the practical interest of 
the previous two approaches, this latter one, referred as fQL, is mainly of theoretical 
interest. The aim in this case is to suggest what might be ultimately captured from 
the quasi-linear theory if we could accurately model the full frequency spectrum of the 
nonlinear saturation, and in particular the portion not captured by the leading normal 
mode (or modes) in each wave-number. The fQL can be interpreted as the limit passage 



from the discrete mode j summation of Eq. (3.37) to an appropriate integration over 
the continuous frequency w, i.e. J^j ~^ I 

As already mentioned, there is no practical way to directly test the full frequency 
spectrum approach, i.e. first finding the quasi-linear weights over a wide range of fre- 
quencies, then secondly capturing the corresponding frequency dependent nonlinear 
spectral intensities for each wave-number. The method here developed to study the 
fQL approach relies instead on the simultaneous treatment of both nonlinear plasma 
species and linear tracers |101| . The simulations retain the main species, i.e. ions and 
electrons at full densities ne,main and ni^main, and identical tracer ion and electron 
species at negligible densities Ue^tr and rn^tr, such that ^J^"'*'] — ^ 1. Hence, 

the main plasma species alone drive the potential fluctuations and the true (not external 
or artificial) saturation spectral intensity. Due to their negligible densities, the tracer 
species have no feedback on the potential fiuctuations through the Poisson equation. 
If both the nonlinear E x B a,s well the linear terms are retained in the gyrokinetic 
equations governing the motion of the tracer species, then the turbulent fluxes of the 
main species and of the tracers are identical. However, if the tracers have only lin- 
ear motion, by artificially deleting the nonlinear terms, the resulting tracer turbulent 

fluxes ^i^^ ,Q{fr^ ,Qi^^ can be truly taken as reliable quasi-linear transport esti- 
mates, here referred as the fQL fluxes. Further details on this method, as it has been 
implemented in the GYRO code, are reported in Appendix [B] The quasi-linear weight 
according to the fQL approach can then be defined only for each toroidal wave-number 



k (while the sum over the index j of Eq. (3.37) has no sense anymore here): 

(3.45) 



QL- weighty 



yfQL i^fQL ^fQL 



The relevance of this method relies on the information carried by the transport ratios 
(overages) between the linear tracers and the nonlinear main species. The overages 

^l%^/'^e^^Q{^r^/Qf^Q{%^/Qe^] s^ould in fact quantify the goodness of the best 
quasi-linear approximation, since no additional hypothesis has been introduced apart 
from the suppression of the nonlinear dynamics in the gyrokinetic equation. It is 
reasonable to expect the overages from the fQL approach to be somewhat smaller 
than those from the mQL approximation. The reason is due to the fact that the fQL 
approach does not discount any less unstable mode, while properly accounting for the 
whole structure of every sub-dominant mode. Conversely, the fQL approach exhibits 
an intrinsic failure of the ambipolarity condition, i.e. r{*^/' ^ r{^^. The tracer species 
dynamics is in fact not self-consistently evolved through the Poisson equation; this point 
is a natural consequence of the nature of the quasi-linear problem, where the feedback 
of the linear particle motions on the fluctuating fields is neglected. Nevertheless, since 
the nearly adiabatic electrons are expected to control the particle transport, the tracer 
fQL ion particle transport is usually neglected. 
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The validity of the quasi-hnear response can then be systematically tested by mean of both 
the mQL and fQL approaches. The key measures of success for the quasi-linear response 
approximation arc then: 

• The quasi-linear over nonlinear overages have to be nearly the same in each transport 
channel (particle, ion and electron energy) and, less importantly, in each wave-number. 

• The overages should be nearly constant across a wide variety of physical case parame- 
ters. 



The validity of the quasi- linear response is tested using the gyrokinetic code GYRO. In 
addition to the G A-standard case already defined in Table |3.1[ two other sets of plasma pa- 
rameters are defined: while keeping fixed all the other parameters, the normalized density and 



temperatures gradients are summarized in Table 3.2 In the G A-standard case (ITG-TEM), 





GA-ITG-TEM 


GA-TEMl 


GA-TEM2 


R/Lti 


9.0 


6.0 


3.0 


R/Ltc 


9.0 


6.0 


3.0 


Rj Ln 


3.0 


6.0 


9.0 



Table 3.2: Density and temperatures gradients lengths characterizing the three sets of plasma 
parameters. 



the ITG mode linearly dominates at low k, while TEM leads for 0.5 < kgps < 0.75. For GA- 
TEMl, one half of the linear modes in the spectrum are ITG and the other half are TEM 
with similar growth rates. Finally in GA-TEM2, TEM linearly dominates at all kgp^ < 0.75. 



Table 3.3 shows the quasi-linear/nonlinear overages for the two-step tests according to the 









jjmQL 1 Jjnl 


GA-ITG-TEM 


17.9/12. 1==1.47 


5.83/3.42=1.70 


-3.21/-2.01=1.79 


GA-TEMl 


21.8/15.0=1.45 


21.0/15.3=1.37 


7.86/5.50=1.42 


GA-TEM2 


43.8/25.9=1.69 


56.3/30.3=1.85 


12.2/5.9=2.07 



Table 3.3: Overages for the two-step tests of the mQL approach as ratios of quasi-linear to 
nonlinear diffusivities for different sets of plasma parameters; diffusivities are in Gyro-Bohm 
units, where xgb = Cs/o-pI- 

mQL approach. The average across the results presented in Table [373] gives an overage of 1.64, 
with an RMS of the normalized deviations of 13.5%, which has to be considered a successful 
result. Fig. |3.6| reports the corresponding overages across wave-numbers as well as the non- 
linear saturation spectral intensity for the three cases. It appears that the overages at higher 
k present, sometimes considerably, deviations. This can be easily understood if recalling the 
cross-phases of Figs. |3.3[ |3.4| and |3.5[ where analogous behavior was found. As previously 
mentioned, this breakdown has little effect on the on the total quasi-linear transport, since 
the most relevant part of fluxes is carried by the low k scales. Recalling the formulation of 
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Figure 3.6: Wave- number spectrum of overage for ion energy (a, e, i), electron energy (b, f, 
j), and particle (c, g, k) transport, and nonlinear spectral intensity (d, h, 1) for ITG-TEM, 
TEMl and TEM2 cases. The horizontal lines indicate the net overage in each channel as 
given in Table |3.3[ 
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Eq. ([3:371) 

(Fig, 



this is due to the weak weight of the nonlinear saturation spectrum at the high k 
3.6). The apparent breakdown of the mQL approach at high k could support the guess 
of a nonlinear wave-particle resonance broadening represented by the additional k^D term 
[311 1102j . In this direction, it is useful to note that the gyrokinetic code GENE has been 
used to show that the E x B nonlinearity can be well represented by a (purely real) linear 
damping in the form kyD for a cold ion TEM turbulence j76| . The really large jump in the 
overages for kePs > 0.5 in the GA-ITG-TEM case can be associate with the transition from 
the dominance of the ion (ITG) to the electron (TEM) branch at higher-fc. The ITG-TEM 
case with its inward (pinched) particle flow, i.e. < 0, seems to have the most tenuous 
validity according to the mQL approach and this point will be deepened in the following. 

Turning to the validation of the full frequency spectrum quasi-linear approach (fQL ap- 
proach) for the GA-ITG-TEM case, it is important to firstly verify the plasma tracers method 
itself. Indeed, when properly retaining the nonlinear terms in the gyrokinetic equations gov- 
erning the dynamics of ion and electron tracers, turbulent diffusivities are found to be iden- 
tical between the main plasma species and the tracer^ 

Table |3.4| provides instead the overages for the fQL approach by deleting the nonlinearity in 
the tracer species. As expected, the tracer particle flow is not ambipolar: the ion tracer parti- 





fQL / nl 


^fQ^ lyril 
A.e,tr 1 A.e,main 




TjfQL 1 Tjnl 

i,tr I i,7nain 


ITG-TEM 


16.7/11.3=1.48 


3.66/3.17=1.15 


-2.7/-1.90=1.43 


-5.1/-1.9 



Table 3.4: Overages for the one-step tracer tests of the fQL approach as ratios of quasi- 



linear to nonlinear; the apparent discrepancy with respect to Table 3.3 for the main species 
diffusivities is due to slightly different grids and time averages. Diffusivities are in Gyro-Bohm 
units. 

cle diffusivity is two times larger than the electron tracer one. Since electron non-adiabaticity 
controls particle flow, we consider only the electron particle fiow overage of 1.43, which is a 
value very close to the ion energy overage. In this case, the electron energy transport presents 
the lowest quasi-linear over nonlinear ratio, with an overage factor which is close to unity. 
As already mentioned, it is reasonable to expect that quasi-linear over nonlinear overages 



obtained according to the fQL (Table 3.4 1 approach are smaller than those from the mQL 
approach (Table 3.3), since the nonlinear spectral intensity is also in part distributed over 
the all sub-dominant, modes which should have smaller quasi-linear weights. 

The validation of the quasi-linear response does not only rely on the structure of the quasi- 
linear over nonlinear overage across the k wave-number and the transport channels. Even 
more importantly for the application to quasi-linear models, these overages should also be rea- 
sonably constant over a wide variation of the plasma parameters. In Figure 3, the GA-ITG- 
TEM case is reconsidered over a wide range of temperatures gradients 6.0 < R/Ln.e < 27.0, 
using both the mQL and the fQL approach. The results of Fig. |3.7| show that the quasi-linear 
over nonlinear overages in both ion and electron energy channels remain nearly constant, even 
though the overall transport level increases by 4-8 times. As expected, the fQL approach 
leads to lower average quasi-linear over nonlinear ratio, about 1.3 versus 1.6 from the mQL 
approach. 



^In particular, for the GA-ITG-TEM case, the results of the GYRO simulation are: Xi,tr/Xi 
12.0/12.0, Xe,tr /Xe,main = 3.22/3.28, De,tr / Dg^rnain = —2.01/ — 2.0. Some Small difference in the electron 
diffusivities are expected, since the tracer electrons (and also both main and tracer ions) are evolved according 
to an explicit numerical method, whereas the main electrons use a mixed explicit-implicit scheme (for fast 
parallel motion) 
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Figure 3.7: Quasi-linear over nonlinear overages across a temperature gradient scan based 
on the GA-ITG-TEM standard case, using both the mQL and the fQL approaches. 



According to the mQL approach, the overage for the particle transport, characterized by 
a strongly inward flow at high gradients (i'e/XGS up to -7.45), considerably breaks away 
from the overages on the energy transport. The mQL particle transport overage fails with an 
underflow at a/Lx — 2.0 (R/Lt — 6.0). This latter breakdown is likely due to the lack of sub- 
dominant modes within the mQL approach. More importantly, the particle transport overage 
shows a monotonic increase with higher temperature gradients, with an over-prediction with 
respect to the nonlinear factor by a factor up to 3. It is worth noting that this failure of the 
quasi-linear particle transport with strong turbulence drives, is consistent with the previous 
findings on the cross-phases S(j) — Sn shown in Fig. |3.4[ where the linear phase shifts signif- 
icantly overestimate the nonlinear ones. The key point that could not be verified through 
the study of the cross-phases is checking if a proper accounting of the linear sub-dominant 
modes can fix this over-prediction and get closer to the nonlinear predictions. The latter 
point can now be investigated thanks to the fQL approach. At high gradients, the quasi- 
linear particle transport over-prediction appears less pronounced (because of the accounting 
of sub-dominant modes), but still present, with quasi- linear over nonlinear ratios 2 times 
bigger than those of the energy transport. Hence, the quasi-linear over-prediction observed 
for high gradients, R/Lt >« 15, has necessarily a different origin, and it can not be related 
to the accounting of sub-dominant modes. This failure is in fact equally evident using both 
the mQL and the fQL approaches, where the latter one correctly retain the whole linear 
modes structures. 

These results are then in favor of recognizing the failure of the quasi-linear response when 
trying to reproduce the nonlinear particle transport of strongly inward flows. At this regard, 
it is worth noting that realistic tokamak plasma conditions foresee rather weak core particle 
source, so that the typical operating point has a tiny or null net particle flow in the core. A 
strong outflow down the density gradient takes place instead at the edge, balancing gas feed 
and wall recycling. The strong inward particle pinch reported in Fig. 3.7 for R/Lt >~ 13, 
is mostly due to ion unstable modes, can then appear as a rather pathological condition. 
For this reason, another temperature gradient scan based on a modified GA-TEM2 case is 
explored. Referring to the parameters reported in Table |3.2[ GA-TEM2 is modified consid- 
ering R/Lti = 3.0 0.3. In this case, the strong density gradient results into a turbulence 
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largely dominated by electron TEM modes while the reduced ion temperature gradient sup- 
presses the ion ITG branch. The particle transport is hence characterized by outward flows 



across the whole scan of the electron temperature gradient 3.0 < R/Ltc < 27.0. Fig. 3.8 
reports the results using the mQL approach for the quasi-linear expectations. Here, the 
quasi- linear/nonlinear overages appear relatively constant for both the ion and electron en- 
ergy but most of all for the particle transport, even at high temperature gradients (compare 
with Fig. 3.7). The results of Fig. 3.8 are then in favor of supporting the guess that the 
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Figure 3.8: Quasi-linear over nonlinear overages across an electron temperature gradient scan 
based on a modified GA-TEM2 standard case, using both the mQL approache. 



quasi-linear response may not appropriate for capturing the nonlinear physics responsible for 
strong inward particle flux, while it appears still reliable for more usual moderate particle 
flows. 



3.3 Summary 

The most relevant findings of the present section are here briefly summarized: 

1. The validation of the hypothesis of the quasi-linear response is systematically and 
quantitatively studied on comprehensive nonlinear and linear gyrokinetic simulations 
of tokamak micro-turbulence. Three different methods are defined and applied: (a) 
characteristic turbulence times, or Kubo-like numbers, (b) cross-phases of the fluctu- 
ating quantities, (c) transport weights. 

2. The Kubo-like numbers evaluated across a wide temperature gradient scan on tokamak 
relevant conditions, fully satisfy the quasi-linear ordering of the characteristic turbu- 
lence times. The effective auto-correlation time is found to be smaller than the eddy 
turnover time, a condition that should ensure that the particles are not trapped in the 
field pattern. 
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3. The phase shifts between the transported quantities and the fluctuating potential high- 
lights a good agreement between the nonlinear and linear (most unstable mode) regimes 
for relevant tokamak parameters. Nevertheless, the cross-phases alone can not be used 
to quantitatively infer the validity of the quasi-linear response. A major caveat concerns 
in fact the accounting of sub-dominant modes. 

4. The properly defined transport weights are finally able to quantitatively address the 
validation of the quasi-linear response. The latter one presents a typical over-prediction 
with respect to the nonlinear transport levels. Most importantly this overall factor is 
shown to be reasonably constant across (a) the different transport channels (particle, ion 
and electron energy) , (b) the variation of the plasma parameters on tokamak relevant 
conditions. Nevertheless, the quasi-linear particle transport appears to significantly 
deviate from the nonlinear expectations when dealing with strong inward flows. 



Chapter 4 



Improving the underlying 
hypotheses of the model 

In the previous chapter, the gyrokinetic quasi-hnear transport modehng has been investigated 
firstly in terms of its formal derivation and actual formulation applied to tokamak plasmas. 
Secondly, the hypothesis of quasi-linear response has been studied and validated by mean of 
systematic comparisons with comprehensive nonlinear gyrokinetic simulations. This chapter 
is instead dedicated to the other fundamental part contributing to the quasi-linear transport, 
i.e. the model of the nonlinear potential saturation, recalling the generic expression of the 
quasi-linear flux Eq. |3.37| Once again we stress the crucial importance of a clear distinction 
between the latter issue and the one examined in the previous chapter. Since the intrinsic 
nature of the quasi-linear theory can not provide any information on the nonlinear satura- 
tion process, the choices and the resulting implications adopted for modeling the saturated 
spectrum deserve careful dedicated analysis. 

The general aim of this chapter is the elaboration of a nonlinear saturation spectrum 
model, which can integrate the information carried by the quasi-linear response in order to 
produce reliable estimates of the turbulent fluxes. 

Following from these preliminary remarks, the tools needed to investigate this issue has to 
be defined. In this work, we propose a novel physics oriented approach, where the information 
is carried in parallel from two different sources: 

1. The properties of the nonlinear saturation process are studied by comprehensive non- 
linear gyrokinetic simulations. As the latter ones represent the most advanced tool 
actually available to gain quantitative information on the tokamak turbulence, the use 
of nonlinear gyrokinetic simulations is essential in order to model the saturation mech- 
anisms. 

2. The experimental fiuctuation measurements are a primal information about the real 
nonlinear processes that are at play in tokamak plasmas. The analysis of these data 
is of great importance in order to validate and improve the hypotheses adopted in the 
model for the saturation spectrum. 

Both the nonlinear simulations and the measurements must equally contribute to the im- 
provement of the potential saturation choices, which otherwise would remain completely 
arbitrary. Conversely, the integration of the experimental and the numerical results allows 



62 



Improving the underlying hypotheses of the model 



relevant advances in the understanding of the underlying physics and therefore justifies the 
choices that will be adopted in QuaLiKiz. 

The approach here just briefiy described, is a unique feature proposed within the present 
work. In fact, in most of the existing quasi-linear tokamak transport models the quasi- linear 
response is coupled to generic mixing length saturation rules, without a clear explanation of 
the underlying choices. More recent quasi-linear transport models, like TGLF, have more 
carefully treated the issue of the nonlinear saturation hypotheses. Nevertheless, the strategy 
adopted by TGLF is different from the one here proposed. In the case of TGLF in fact, 
the quasi-linear gyro-fluid response is coupled to an artificial nonlinear saturation intensitjj^ 
which proves superior when fitting the resulting TGLF total fluxes with a large database of 
GYRO nonlinear simulations. 

Conversely, in the present work we do not rely on a numerical database best flt approach, but 
the fluctuation measurements are used to physically interpret and validate the expectations 
provided by the nonlinear gyrokinetic simulations. Once we trust these first principle numer- 
ical results, they can be used to powerfully improve the choices retained in the saturation 
model of a quasi-linear transport model. 

The chapter is organized as follows. The first part addresses the issue of the validation of 
the numerical tools, i.e. the nonlinear gyrokinetic simulations, against the turbulence mea- 
surements done in the Tore Supra tokamak. More in particular, after a brief introduction 
about the experimental techniques used to diagnose the tokamak plasma fluctuations, de- 
tailed analysis is devoted to the spectral structure of the density fluctuation in terms of both 
the wave-number k and the frequency a; dependences. Consequently, the choices concerning 
the saturated potential adopted in QuaLiKiz are discussed. Here, both the spectral structure 
(in the k and co spaces) and the saturation rules adopted to weight the different contributions 
to the total transport are treated. 



4.1 Validating the nonlinear predictions against mea- 
surements 

The understanding of turbulent fluctuations in the core of tokamak plasmas, causing a degra- 
dation of the confinement, is a crucial issue in view of future fusion reactors. Global empirical 
scaling laws based on the existing experiments [llOi [55] are often used to extrapolate the 
performance of the next generation devices like ITER and DEMO. However, predictive ca- 
pabilities should rely on first principle models retaining comprehensive physics. Hence, the 
development of validated predictive codes is an essential task for ensuring the success of the 
present and the future fusion experiments. 

While the general tokamak plasma instabilities framework has an history of intensive ana- 
lytical and computational efforts from more than 30 years, the recent times are characterized 
by the following very relevant advances: 

• Numerical codes matured to the point to realistically include what it is believed to be 
the essential complexity of tokamak physics, allowing to provide quantitative informa- 
tion about the plasma turbulence and the associated transport. 

• Great progresses have been achieved in developing reliable diagnostics, able to accu- 
rately measure the properties of the tokamak turbulence in wide range of parameters 
and experimental conditions. 



^called V-norm, further details on this point can be found in I59| 
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• The computing power has significantly raised, especially thanks to the massive use of 
parallel high performance computing, allowing to perform highly sophisticated simula- 
tions on tractable time scales. 

The issue of the nonlinear codes validation has very recently produced a great amount of re- 
search. Some of the most relevant efforts have been done across several tokamaks worldwide, 
namely on DIII-D [SaHUH], Alcator C-mod [53 [SH] and Tore Supra [23]. 

In order to truly validate a tokamak turbulence model, two different levels of information 
have to be accurately predicted. The first one refers to the detailed spectral structure of 
the turbulence (lower-order quantities), i.e. both the spatial and time spectral shapes and 
amplitudes of the different fluctuation fields; the second level instead deals with the macro- 
scopic observables (higher-order quantities) , such as the total turbulent flows and fluctuation 
amplitudes, that integrate the previous spectral quantities. This approach has also been 
summarized by a practical strategy towards the validation in the fusion research, identifying 
[95] the following two key concepts: 

1. Primacy Hierarchy: Ranking of a measurable quantity in terms of the extent to which 
other effects integrate to set the value of the quantity. Assesses ability of measurement 
to discriminate between different non-validated models. 

2. Validation Metric: A formula for objectively quantifying a comparison between a sim- 
ulation result and experimental data. The metric may take into account errors and 
uncertainties in both sets of data as well as other factors such as the primacy of the 
quantities being compared. 

In this section, nonlinear gyrokinetic simulations are quantitatively compared to turbulence 
measurements on the Tore Supra tokamak. The relevance of this analysis relies on the fact 
that the high-order scalar observables are coherently verified also through the investigation 
of the lower-order spectral quantities. With respect to the previous works, the present study 
is qualified by the simultaneous validation of (1) the total heat transport coefficient, (2) rms 
values of the density fluctuations (5n, (3) kg and (4) kr wave-number Sn spectra and (5) the 
frequency Sn spectra. 

4.1.1 Turbulent fluctuation measurements 

Several techniques have been developed in order to investigate the spectral structure and 
the main features of the density and temperature fluctuations in tokamak plasmas. An ex- 
haustive review of these diagnostic systems is largely beyond the scope of this work. In Fig. 
|4.1[ we just remind some of the actual most important techniques, comparing their general 
spectral range in k± to the typical spatial scales of the tokamak instabilities. The diagnostics 
appearing in Fig. |4.1| are (details can be found in the cited references) : the Beam Emission 
Spectroscopy BES |75). the reflectometry (Doppler [3S], fast-sweeping and fixed frequency 
[^[39]), the Far Infrared FIR [16] and laser scatterings [SZlllI], the Correlation Electron 
Cyclotron Emission CECE [108]. 

The TORE SUPRA tokamak is particularly well suited to study the local spectral struc- 
ture of the density fluctuations in the low and medium range k±ps < 2: it is in fact equipped 
with complementary microwave diagnostics, fast-sweeping [251 and Doppler [48j reflectome- 
ters. Though the geometrical arrangement and measurement techniques are different, these 
two methods are both based on the detection of the held backscattered on density fluctua- 
tions, whose wave number matches the Bragg rule kf = —2ki, where ki is the local probing 
wave-vector. Synthetic schemes of both the reflectometry systems are presented in Fig. |4.2| 
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Figure 4.1: Picture summarizing some actual tokamak fluctuation measurement techniques, 
comparing the diagnostic typical spectral ranges to those of the ITG, TEM and ETG insta- 
bilities. 



In the fast-sweeping system, the probing wave is launched in the equatorial plane, with 
a ki in radial direction. The fast sweeping of the probing frequency (from 50 to 110 GHz 
in 20/is) allows sensitive measurements of the density profile n{r). This acquisition time is 
sufficiently fast and at high repetition rate, in order to extract the density fluctuation profile 
5n{r). The latter is obtained from the integration of local spectra S[Sn/n\ {kj.,r) between 
1 < fcr < 10 cm""'^. The primal quantities measured by the refiectometer are the phase 
fluctuations, which are linked to the plasma density fluctuation through a transfer function: 
an iterative process of comparisons between full- wave simulations and experimental measure- 
ment is used (for more details see Ref. [39]). The local kr spectrum can then be obtained 
using a sliding spatial Fourier transform of 6n/n. The typical size of the spatial window 
allows probing the spectrum between 1 ^ A:^ ^ 20 cm~^, with a resolution on the spatial 
localization of Ar/a ~ 0.04 (a is the plasma minor radius). Note that in this method, fluctu- 
ations with poloidal wave- numbers up to kg « 10 cm~^ are also included in the signal due to 
flnite beam size. This has to be taken into account when comparing with theory or simula- 
tions. The probability density function of Sn/n (r, t) shows a log- normal distribution law for 
all radial positions. Then time averaged 6n/n (r) values are obtained over a statistics of 1000 
measurements. The computation of transfer function uses a process based on 100 normally 
distributed random samples, implying thus a 10% uncertainty on its determination. The 
uncertainties coming from the sensitivity to the temperature and density are also accounted. 
Finally, the uncertainties on Sn/n (r) are estimated by quadrature addiction, with a typical 
level of 15 %. 

The Doppler reflectometry technique shares with the standard reflectometry the advan- 
tages of good spatial and temporal resolution, easy access and low cost. Differently from 
the fast-sweeping reflectometry, in the case of the Doppler reflectometry system, the prob- 
ing beam is launched in oblique incidence with respect to the cut-off layer. The scattering 
process is mostly localised at the cut-off, whose position is set by the (fixed) frequency of 
the beam |3Sj. In the case of 0-mode beam polarisation, the selected wavenumber at the 



4.1 Validating the nonlinear predictions against measurements 



65 



BACKSCATTERIKG DT Sn/n 



ko 




BACKSCffTTCRING 

BYHIGHkran/n 



^ BACKSCATTERED 
SIGNALS 

WAVENUMBER 




cut-off layer 
~ few cm 



radial position 



b) 



Figure 4.2: Schematic representation of (a) the perpendicular incidence reflectometry and 
(b) the Doppler reflectometry. 



cut-off is mainly poloidal. The rms value of the signal is directly proportional to the power 
spectrum of density fluctuations at the selected wave-number: the wave-number spectrum is 
then obtained by varying the antenna tilt angle. This allows us to probe the spatial domain 
0-5 ^ "r/a < 0.95 and the range of wave-numbers 4 < /eg < 15 cm^^, where only very small 
kj. < 1 cm~^ are included in the signal. 



Experimental set-up 



One Tore Supra L-mode ohmic discharge, TS39596, has been chosen as target for the 
validation of the nonlinear simulations against the turbulence measurements. TS39596 is a 
Tore Supra representative standard shot, characterized by high reliability and repeatability. 
This simple discharge in fact does not differ from many other ohmic shots and it presents a 
large time interval with steady plasma profiles and no external momentum input. The time 



evolution of the principal plasma parameters is shown in Fig. 4.3 The plasma parameters 




Figure 4.3: Time evolution of the plasma current, central electron temperature, line averaged 
density and diamagnetic energy content for the discharge TS39596. 
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Table 4.1: Main plasma parameters for the discharge TS39596; electron and ion temperatures 
are referred to r/a — 0.5 



characterizing TS39596 are moreover summarized in Table 4.1 



The set of diagnostics installed on Tore Supra and here used to measure the plasma 
quantities is here briefly described. The radial profiles of the electron density are measured 
by the fast-sweeping reflectometers [53] for magnetic fields Bt > 2.82 T, since the localization 
of the reflectometry cut-off layer, and hence the measurements, depends on the values of the 
plasma Bt- Since TS39596 has Bt — 2 AO T, the electron density profile is obtained by an 
Abel inversion of the line density measurements made by the interferometry. The resolution 
on the radial location of the measurements is 5 cm when using the interferometry chords, 
and 1 cm when using the reflectometry. The uncertainty on the density values is estimated 
at 10%. 

The Radial electron temperature profiles are measured by both Thomson scattering and 
electron cyclotronic emission ECE [89]. The radial location resolution of the measurements 
is 10 cm for the Thomson scattering, it is while close to 1 cm when using the ECE. The 
measurements uncertainties are estimated at around 15% for the Thomson scattering and 
7% for the ECE. 

The ion temperatures are measured by the charge exchange diagnostic |49) . Eight chords are 
available up to normalized radii of r/a = 0.7. The experimental uncertainty is around 20% 
for the core region and 15% in the gradients zone. The radial resolution varies between 2 
and 6 cm. High quality temperature and density profiles are crucial in order to minimize the 
resulting uncertainties on the gradient lengths, which are derived from radial derivatives of 
these quantities. 

The global effective ion charge Z^f / is measured by Brehmstralhung emission. The estimate 
on the radiative power given through the Matthews law by the line integrated measurements, 
is consistently checked with the direct radiative power measurements. The 2^e// radial profile 
can not be obtained in this case, thus a flat Z^fj profile is assumed. Simulations done with 
the CRONOS code [3] using the measured Z^ff have verified that the current diffusion is 
compatible with the measured flux consumption. 

Finally, the safety factor q (r) profile is also estimated using the CRONOS runs constrained 
by 10 polarimetry angles measurements. The current diffusion reproduces both the measured 
flux consumption and the internal value of inductance. 

The discharge TS39596 has been analyzed by mean of nonlinear gyrokinetic simulations 
on 4 different radial locations, namely r/a = 0.4, r/a = 0.5, r/a = 0.6 and r/a = 0.7. The 
experimental parameters characterizing the plasma at each location are reported in Table 
4.2[ The dimensionless parameters appearing in the Table [T2| are particularly worth of note. 
These are p*, the ratio of the the ion-sound Larmor radius and the plasma minor radius, /3, 
the ratio of the plasma kinetic and magnetic pressure, and , the electron- ion collision fre- 
quency normalized to the electron bounce frequency. The last ones are generally recognized 
as relevant parameters for fusion plasmas, so that global scalings of the energy confinement 
time have been proposed with crucial dependence on these three dimensionless parameters. 
Despite the low value of the magnetic field, the discharge TS39596 exhibit p^, dimensionless 
parameters which are not so far from the values expected for ITER. At the same time, the 
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Table 4.2: Experimental parameters for the discharge TS39596 at different radial locations; 
p* = Ps/a, i^ei frequencies are in units of c^/a (c^ is the ion-sound speed) and s is the magnetic 
shear. 

collisionality in TS39596 is very high, while in ITER very low i^^ values are foreseen. In fact, 
adopting the following definition, i^^, — v^i ^/^^/t i where e = t/Rqis the local inverse aspect 
ratio, for TS39596 we obtain = 0.66, 0.82, 1.05, 1.60 respectively for r/a = 0.4, 0.5, 0.6, 0.7. 

Local nonlinear gyrokinetic simulations 

Local nonlinear gyrokinetic simulations on TS39596 are performed using the GYRO code. 
The latter one has in fact the capability to treat several among the key physical mechanisms 
needed to gain a quantitatively relevant information on the tokamak turbulence. These 
nonlinear simulations are performed in the local (flux-tube) limit of vanishing p* , considering 
the 4 radial locations r/a — 0.4, 0.5, 0.6, 0.7. The GYRO simulations make use of all the 



experimental parameters reported in Table 4.2 without any modification. A summary of the 
numerical grid adopted is here given. 

• Electromagnetic fluctuations are retained. 

• Collisions are treated trough the electron-ion pitch-angle scattering operator. 

• Realistic magnetic geometry via the Miller equilibrium is adopted (i.e. flnite aspect 
ratio effects are considered). 

• Gyro-kinetic ions and drift-kinetic electrons are retained (i.e. electron FLR effects are 
neglected) . 

• Zf,ff = 1, hence rig = rii considering the real mass ratio yjra^ml — 60. 

• Box size in the perpendicular directions [Lx/ Ps, Ly / ps] = [122, 122] with a radial reso- 
lution Ax/ps = 0.30 (i.e. 400 radial grid points). 

• 128-point velocity space discretization per spatial cell: 8 pitch angles, 8 energies and 2 
signs of velocity. 

• 12 grid points in the parallel direction, 18 in the gyroaverage and 5 in the radial 
derivative. 

• The simulation at r/a — 0.7 use 32 complex toroidal harmonics resolving 0.0 < kgps < 
1.59 scales, while for r/a = 0.4, 0.5, 0.6, the range 0.0 < kgps < 1.0 is solved with 20 
Fourier modes. 
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• Statistical averages on the saturated nonlinear state are performed on a typical physical 
equivalent time of 2.6 ms. 

Nonlinear GYRO simulations do achieve statistically steady states of turbulence through 
dissipative upwind advection schemes in the real (not velocity) space. All the simulations 
presented in this thesis work verify the condition that both the flux and the fluctuation levels 
are well converged with respect to the velocity-space resolution at fixed spatial resolution. 
The velocity space resolution here employed (128 points for spatial grid point) is the standard 
one for typical GYRO simulations: the convergence with respect to this standard velocity- 
space resolution is systematically investigated in [201. This aims to check that there is no 
missing velocity-space structure that would affect either transport levels or entropy produc- 
tion. Indeed this kind of approach is based on the statement that spatial upwind numerical 
schemes, hence finite numerical dissipation, provide a physically meaningful procedure to 
achieve turbulent steady states. Obviously the strategy followed by the GYRO code is not 
the only one possibility; other solutions deal with the issue of the dissipation looking for the 
direct implementation of a physical collision operator, like the recent proposals by [TJ[S]. 

It is important to briefiy discuss the two main physical limitations of this kind of simu- 
lations, namely the local approach and the lack of electron spatial scales. This first issue is 
not expected to be a severe limit. Even if the argument is still an open subject of research, 
it is widely recognized that for such small values of w 2 • 10"'^, the local treatment of the 
turbulence appears as a very good approximation. This has been also demonstrated by the 
success in recovering the breaking from the gyro-Bohm (corresponding to the local limit of 
small p*) to the Bohm transport scaling using nonlinear gyrokinetic simulations |19j . On 
the other hand, the fact that the electron spatial scales, i.e. kgpg ^1.0 are here neglected, 
deserves some additional care. The main reason for this simplification is the extremely high 
computational cost required by the coherent treatment of this multi-scale nonlinear problem. 
The problem is still presently an open issue. A first example of massive nonlinear gyrokinetic 
simulations of coupled ITG-TEM-ETG turbulence |100j indicated small relative contribu- 
tions to the total transport level coming from the high-A: scales. More recently instead, in 
Ref. [H], it has been shown that in the case of reduced low-fc linear drives, a scale separation 
between electron and ion thermal transport is observed, originating a non negligible trans- 
port contribution at the small electron scales. Nevertheless, both these example do not adopt 
the real ion/electron mass ratio, because of computational cost constraints. Therefore, the 
present analysis focuses only on the spatial scales smaller or comparable to the ion Larmor 
radius, i.e. kgpg < 1.0, where the low k ranges are well resolved by the nonlinear simulations 
and the spectral quantities provide a quantitatively relevant information. 

Results on the higher-order quantities 

The primal experimental quantities that have to be compared to the expectations by the 
nonlinear gyrokinetic simulations are the higher-order scalar quantities. Here, we focus on 
both the total heat transport coefficients and the RMS values of the density fluctuations. 
The total effective heat diffusivity Xeff is experimentally obtained from a power balance 
analysis, performed with an interpretative CRONOS run, with the following definition: 

where qe and qt are the electron and ion heat fluxes respectively. 

The experimental uncertainty is estimated taking into account the time evolution of the pro- 
flies during 1 s. Unfortunately, it is not possible to have a quantitatively reliable information 
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on the separate contribution from the electron and ion heat transport; in the case of high 
coUisional plasmas like TS39596 in fact, that would imply a detailed knowledge on the cou- 
pling terms between the two species. 

It is worth noting that in the case of steady plasma conditions, as verified for TS39596, the 
plasma profiles accommodate to null particle flows in the core, because no source is present 
in that region. For this reason, in absence of convection terms, the heat flux calculated by 
the CRONOS interpretative run effectively coincides with the energy flux. Conversely, the 
gyrokinetic simulations do not solve the physical problem at fixed fiux, whereas the temper- 
ature and density gradients are the quantities that are imposed. A certain particle flux is 
then calculated by the code and has to be appropriately treated. At this scope, the following 
relation is used, that links the energy and the heat fluxes for a generic species s: 



Qs 



:r,r.. 



(4.2) 



In Eq. (4.2), Qg and Tg are respectively the energy and particle fluxes predicted by GYRO, 



averaged over a flux-surface and a nonlinear saturation time interval (see Appendix |B] for 
details). 

The experimental results are compared with the expectations from the GYRO simulations 



in Fig. 4.4 A very good agreement between the numerical expectations and the experimental 
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Figure 4.4: Radial profile of the experimental effective heat diffusivity and comparison with 
the GYRO predictions. 



profile of Xef f within the error bars is achieved. The exception of the radial point correspond- 
ing to r/a — 0.4 is due to marginal turbulence found by the simulation on the experimental 
parameters. Nevertheless, the stiffness to V^Ti e (whose experimental uncertainty is about 
±30%) inherent in the transport problem |19| . may induce to find an agreement with the 
experimental findings if varying the normalized gradients within the uncertainty levels. This 
argument has in fact been demonstrated already in several studies [1^1 |Ml [M] ■ This kind 
of considerations should suggests that a reliable validation of the turbulence model can not 
be limited to the comparison of the scalar Xef / quantity. The approach here adopted relies 
instead on a different strategy for the validation of the nonlinear simulations, which are con- 
strainted to the simoultaneous comparison with several scalar and spectral measurements. 
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The radial profile of the RMS density fluctuations is provided by the fast-sweeping re- 
flectometry and is shown in Fig. |4.5| In order to validate the predictions by the nonlinear 
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Figure 4.5: Radial profile of the experimental RMS 5n/n and comparison with the GYRO 
predictions. 



simulations, it is essential to properly reconstruct the same quantity that is measured by the 
diagnostics. Recalling the previous description of the fast-sweeping reflectometry system, the 
following relation is applied to the GYRO results: 
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(4.3) 



In Eq. (4.3), Sn/ii {kg, k^) refers to the physical density fluctuations predicted by the GYRO 
simulations computed at the outboard midplane, according to the diagnostic's line of sight. 
For simplicity of notation, the time average over the nonlinear saturation phase is skipped 



in the equation. The results shown in Fig. 4.5 reports a remarkable quantitative agreement 
with the experimental radial 5n/n profile within the error bars. Also the slight increase of 
the density fluctuations towards the external radii found by the diagnostic is matched by the 
GYRO simulations. 



4.1.2 The wave- number spectrum 
Results on the lower-order quantities 

The local nonlinear GYRO simulations have been shown to simultaneously reproduce the 
radial profllcs of both the total effective heat diffusivity Xef f and the RMS density fluctua- 
tions Sn/n. In order to further verify the effectiveness of the turbulence predictions by the 
nonlinear gyrokinetic code, it is important to examine the spectral structure of the fluctua- 
tions, possibly in both the perpendicular directions to the magnetic field. 

As already discussed, the large amount of free energy in tokamak plasmas originates a 
wide variety of unstable modes. Due to the intense magnetic field, those are anisotropic: 
the convective cells exhibit much smaller parallel (to the magnetic field) wave vectors than 
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transverse ones, fc|| <^ fcj_, such that tokamak turbulence is quasi two-dimensional [35]. The 
amount of turbulence anisotropy in the perpendicular plane is a particularly relevant issue for 
tokamak applications: any asymmetry favoring the formation of radially elongated structures 
could in fact enhance the cross- field anomalous transport, with consequent loss of confine- 
ment. Analysis of the fluctuation spectral power density in the transverse wave-number space 
{kg,kr) (poloidal and radial wave-vector, respectively) allows to characterize the turbulence 
structure and gives insight into its dynamics in terms of flow of the turbulent energy E{k) at 
the scale l/fc. Following Kraichnan [51], the injected energy in two-dimensional fluid turbu- 
lence is expected to be nonlinearly distributed through an inertial range, leading to a direct 
enstrophy cascade E{k) — k~^ down to the dissipative length, and an inverse cascade towards 
larges scales, with E{k) = k~^^^. Nevertheless in tokamak turbulence, three-dimensional ef- 
fects, multi-field dynamics and energy injection at disparate spatial scales may originate 
significant departures from this simplified picture, thus requiring dedicated nonlinear simu- 
lations. On most experiments (SSI [1^1 HT] the density fluctuation wave-number spectrum 
S{k±) = \Sn{k±)/n\'^ shows a decay S{k±) ~ /c" with a = —3.5 ± 0.5, whilst assuming an 
adiabatic electron response, Kraichnan dual cascade would predict a — —14/3 and a = —6, 
respectively for the inverse and forward ranges. Kraichnan 2D dual cascade foresees in fact 
the energy density scalings E (fc) ^ k^^^^ and E (fc) ~ fc"-^. In 2D, the energy in an infinites- 
imal shell of thickness dk can be written as E (k) ^ J {■^mv'^^ kdk. Since Vk ^ i '"^'' , one 
finds E (k) ^ k^ \d(j)k\^. The scalings for the fluctuating potential are then |(5(/)fc|^ ^ k~^'^^^ 
and ^ k~^. These scalings are only here presented to give a reference that applies 

to the more common 2D fluid turbulence, which anyway cannot be applied to the case of 
tokamak plasmas. 



Among the 4 different radial positions considered by the GYRO simulations (see Figs. 4.4 
and 4.5), here we focus on the local nonlinear simulation at r/a ^ 0.7. The latter is used to 
validate the gyrokinetic predictions against the measurements of the bidimensional {kg,kr) 
density fluctuation spectra. Because of the low value of the magnetic field in the discharge 
TS39596 in fact, the data from the fast-sweeping and from the Doppler reflectometry are 
simultaneously available only for r/a > 0.7. 

The density fluctuations kg spectrum at r/a = 0.7 from the Doppler reflectometry is pre- 
sented in Fig. |4.6[ Focusing on the region kgps < 1.0, the experimental data exhibit a power 
law decay with spectral index ag = —4.3 ± 0.7. Possible transition towards steeper slopes at 
kgps > 1 as reported in Ref. |47j is not addressed here. In order to correctly reproduce the 
quantiy measured by the Doppler measurements, the following relation is adopted: 
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which is used on the density fluctuations simulated by GYRO at the outboard midplane, ac- 
cording to the diagnostic's line of sight. The simulation gives a spectral index of ag = —4.3 
for 0.4 < kgps < 1.0, in remarkable agreement with the reflectometry data within experi- 
mental uncertainty. 

An additional insight on the kg spectrum of the turbulent density fluctuations has been 
gained also by mean of global gyrokinetic simulations performed with the semi-Lagrangian 
code GYSELA. The latter one retains in fact the radial variations of the background gradi- 
ents, conversely to the GYRO local simulations reported here above. Since the equilibrium 
profile is self-consistently evolved in such a code, the equilibrium gradients decay with time 
due to turbulent transport. Long enough simulation runs allow the system to reach a new 
equilibrium, characterized by well defined averaged profiles which are controlled by both the 
prescribed temperatures at the radial boundaries and the turbulent transport level. The main 
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Figure 4.6: a) Experimental density fluctuation kg spectrum at r/a — 0.7 from the Doppler 
reflectometry (squares), and comparison with the GYRO (diamonds) and GYSELA (trian- 
gles) predictions; b) Comparison between the Doppler kg spectrum and GYRO simulations 
retaining different physics, namely with Miller geometry - kinetic electrons - electromag- 
netic fluctuations {full physics, diamonds), with circular s — a geometry - kinetic electrons 
- electrostatic fluctuations (circles) and with circular s — a geometry - adiabatic electrons - 
electrostatic fluctuations (down-triangles). 
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additional difference with respect to the GYRO simulations is the adiabatic assumption for 
the electron response, such that density and electrostatic potential fluctuations are equal, 
i.e. dn/n = ebc^jT,^. This approximation can be though justified in this case due to the high 
levels of coUisionality, with a consequent reduction of the non-adiabatic electron response 
coming from the trapped electron modes. The local Tore Supra parameters of Table [42] have 
been matched in the global simulations at rja — 0.7, but the normalized gyroradius was 
increased up to p* = 8 • 10"'^ because of limited numerical resources. It should be noticed, 
however, that such a mismatch is likely not to impact the results provided the turbulence 
exhibits a gyroBohm scaling, as expected at these low values. 

A full torus simulation on 0.5 < rja < 0.9 has been run with GYSELA using the grid points 
resolution [r, 9, ip, wy , fj] = [256, 256, 128, 64, 8] {9 and ip are the poloidal and toroidal angles, 
W|| the parallel velocity and ^ the adiabatic invariant). GYRO makes use of a field-aligned 
coordinate system, while GYSELA operates with the poloidal 9 and toroidal ip angles. For 
this reason, specific relations on the different Fourier expansions of the two codes, have to be 
adopted in order to compute the same spectra at = needed for the comparison with the 
measurements. These details are reported in Appendix [B| 

The GYSELA kg fluctuation spectrum is also shown in Fig. |4.6[ giving a spectral expo- 
nent ag = —5.2 for 0.4 < kgps < 1.0. The purpose for which the GYSELA results are here 
shown is not the same as for GYRO. The aim of the GYRO simulations was to demonstrate 
that several measurements can be quantitatively and simultaneously matched when retain- 
ing as much physics as possible in the underlying model: the two main limits in the case 
of GYRO are the local approach (justified by the low « 2 • 10"'^) and the cutoff on the 
small electron scales (whose impact on the fluxes should be sub-dominant and moreover these 
small electron scales are not resolved by the reflectometers in this discharge). The actual 
capabilities of GYSELA do not allow such a multi-level quantitative comparison with the 
experimental data. The reason why the GYSELA |(5(/)fce|^ spectrum is here shown follows 
from the hypothesis that the tokamak turbulence wave-number spectrum is dominated by 
the E X B nonlinear convection terms, providing a rather general shape for the kg spectrum. 
Even if the similarity of the GYRO and GYSELA spectra for kgps > 0.4 can go in this direc- 
tion, of course this statement still remains an hypothesis. The difference between the GYRO 
and GYSELA spectra at low k is mostly due to the not small differences in the underlying 
model. Due to the intrinsic global nature of GYSELA (hence numerical cost), the p* used 
in the simulation has been increased by a factor of 5 with respect to the experimental value 
(p* = 1.6 — >■ 8 • 10~^): large sheared flows in the range of the geodesic acoustic modes are 
damping the large scale fluctuations, hence flattening the spectrum. For this reason and the 
lack of the electron physics in the model, the transport resulting from the GYSELA simula- 
tion cannot be compared to the expectations by GYRO. 

The latter argument could rise the following question: which is the essential physics that 
has to be retained in the gyrokinetic simulations in order to match the measurements within 
the experimental uncertainty? This point is investigated again by mean of GYRO simula- 
tions, retaining different levels of complexity in the underlying physical model. Fig. |4.6[ b) 
reports the results. The first GYRO case deals with the full physics simulation, i.e. using 
a proper Miller magnetic equilibrium coupled to the treatment of kinetic electrons and elec- 
tromagnetic fluctuations. At lower complexity, two other simulations are performed, both 
adopting the simplified s — a geometry: hence a second one still includes the kinetic electrons, 
while the last one retains only adiabatic electrons. It is worth noting that in principle, this 
latter case should correspond to the local (p* — >■ 0) limit of the global GYSELA simulation. 
The results shown in Fig. |4.6[ b) highlight that, even if differences are present, the three 
GYRO spectra retaining different physics, still match the Doppler measurements within the 
experimental uncertainty in the range 0.4 < kgps < 1.0. The power spectral indexes varies 
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from the ag = —4.3 for the full physics simulation to ag = —5.1 in the case of the s — a 
adiabatic electron case: as already seen by GYSELA, this confirms the expected result that 
neglecting the non-adiabatic electron response produces steeper spectral slopes. Moreover, 
the adiabatic electrons simulations (both GYRO and GYSELA) lead to an apparent better 
agreement on the high-fc part of the spectrum; nevertheless this effect is only a fortuitous 
effect. It appears that the most relevant differences when changing the physical complexity 
in the GYRO simulations, are present at the low-A; scales kgps < 0.4. These changes have a 
non negligible impact on the total transport level, which is in fact mostly carried by the scales 
kgps ~ 0.2. The simplified simulations evidently underestimate the total level of turbulence 
with respect to the full physics case, which has instead already demonstrated to correctly 



match the experimental values of the heat flux (Fig. 4.4 1. On the other hand, the turbulence 
diagnostic has not access to the kgps < 0.4 scales; hence in this case, the Doppler reflectom- 
etry alone can not be used to discriminate between the different physical effects retained in 
the simulations. 

The local kr density fluctuation spectrum is obtained by the fast-sweeping reflectometry; 



the measurements at the same radial position r/a = 0.7 for TS39596 are shown in Fig. 4.7 
The measurements still exhibit a power law decay with a spectral exponent = —2.7 ±0.25 
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Figure 4.7: Experimental density fluctuation kr spectrum from the fast-sweeping reflectom- 
etry at r/a = 0.7 , and comparison with the GYRO predictions. 



for scales corresponding to 0.4 < krPs < 2.0. This spectral quantity is reconstructed through 
the relation: 
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(4.5) 



which is applied to the GYRO Sn predictions at the outboard midplane. A very good 
agreement with the fast-sweeping reflectometry data is achieved, both in the magnitude 
and the slope of the relative fluctuation level, covering also the larger spatial scales up to 



^r,mm ~ 1 Cm 
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Discussion on the 6n k spectra 



Recalling the results on the wave-number density spectra just presented in Figs. |4.6| and 
4.7 it clearly appears that at the same radial location r/a = 0.7, the two reflectometers 
provide different fluctuation spectral exponents in the perpendicular plane. This discrepancy 
is above the experimental uncertainties, showing in particular ag = — 4.3±0.7 by the Doppler 
reflectometry while = — 2.7±0.25 according to the fast-sweeping reflectometry, both in the 
range of k±ps > 0.4. Such a discrepancy could initially suggest a highly anisotropic turbu- 
lence, favoring the formation of radially elongated structures. The nonlinear simulations can 
be used in this case as a powerful tool for clarifying this relevant issue. The GYRO results 
motivate in fact a revised interpretation of the experimental evidence: the two dissimilar 
exponents may be simply ascribed to intrinsic instrumental effects. 

In the case of the fast-sweeping kr spectrum, the contributions of medium-low kg < 
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Figure 4.8: (a) \6n/n\'^ spectra from the GYRO simulation at r/a = 0.7, and the impact 
of reconstructing the Doppler reflectometry instrumental response on the kg spectrum, (b) 
Contour plot of the GYRO logiQ\5n{kr, kg)/n\'^. 
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10 cm~^ wave-numbers are retained, as it appears from Eq. (4.5). On the other hand, 
with the Doppler reflectometry, the kg spectrum selects only very low radial wave-numbers 
kr < 1 cm~^, as exphcated by Eq. (4.4) [48]. The density fluctuation spectra predicted by 
the GYRO simulations reflect this asymmetry. More in particular, as it is shown in Fig. 
|4.8[ the kg spectral exponents given by GYRO clearly exhibit a difference when integrating 
over the diagnostic Doppler range < fc^ < 1 cm"-'^ rather than accounting for all the radial 
wave-numbers. In the first case in fact, ag — —4.3 is obtained, the results already shown in 
Fig. |4.6| and in agreement with the measurements. In the second case, the spectral ag w —2.9 
is produced. 

The numerical predictions observe that the strong anisotropy carried by the peak in the kg 
axis significantly affects the k±ps < 0.4 ranges, while the asymmetry appears weaker, but 
still present, at smaller spatial scales |1001 [55] , The iso-level contours of \Sn/n{kj., kg)]"^ in 
the perpendicular plane computed by GYRO and reported in Fig. 4.8 b) confirm this kind of 
expected picture, identifying a linearly driven turbulence anisotropy around the k±ps ~ 0.2 
scales, which are not accessible by our Doppler reflectometry system. 



4.1.3 The frequency spectrum 

Firstly, it is of interest to compare some quantities concerning the frequency of the turbulent 
modes provided by the numerical simulation on the nonlinear and the linear regimes. The 
fcg-dependent frequency spectrum computed by the nonlinear gyrokinetic simulation pro- 
vides a relevant information on the dominant nonlinear frequency w^'' and broadening Aw^?' 
at a given wave-number. Both these quantities can be compared respectively to the linear 
frequency a;[™ and growth rate 7^™. This kind of study is directly useful to validate and 
impro ve the choices assumed in the quasi-linear transport models. As discussed in the para- 
graph 3.1.2 QuaLiKiz keeps the linear frequencies in the calculation of the turbulent 



fluxes, while it assumes a frequency broadening with the widths equal to the linear growth 
rates, Awfc = 7^. 

The results reported in Fig. |4.9| show that for kgps < 0.5, where most part of the 
transport is driven, the mean nonlinear frequencief]^ are in good agreement with the linear 
ones of the most unstable mode. For kgps > 0.5, the linearly most unstable mode jumps to 
the electron diamagnetic direction with increasing while nonlinearly o;^?' tends to zero, 

with a non negligible dispersion. On the other hand, the nonlinear frequency broadenings 
show a departure from the growth rates of the linear most unstable mode. The numerical re- 
sults illustrate that the nonlinear decorrelations operate stronger than the linear ones, hence 



Aw]!' > 7^*". The evidence of Fig. 4.9 is then consistent with the previous evaluation of 



the Kubo-like numbers from nonlinear gyrokinetic simulations (in paragraph 3.2.1 Fig. 3.2). 
This point certainly deserves more accurate study in the future; interestingly, even if this 
observation is at odd with the actual choices of QuaLiKiz, this seems not to dramatically 
affect the total turbulent fluxes, which will be shown in the following of this work to correctly 
reproduce the nonlinear expectations in a number of cases. 

From the experimental point of view, the investigation of the frequency spectra of the 
density fluctuations, is not a very recent subject of research. Relevant works documented 
in the literature are for example Refs. [73l El [48]. Most part of the theoretical efforts on 
this subject are related to the statistical theories of turbulence, like the DIA [1021 IM] , On 



^The moan nonlinear frequencies are computed from the simulation as the ones corresponding to the 
maximal amplitude of the frequency spectrum at a given wave- number kg; each value has an associated 
statistical variance. 
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Figure 4.9: (a) ke spectrum of the frequencies of the linear most unstable mode and of the 
nonlinear frequencies (the bars indicate the statistical variance from the mean value), (b) 
kg spectrum of the linear growth rates and of the nonlinear frequency widths. The results 
refer to nonlinear and linear (most unstable mode) simulations usuing GYRO on TS39596 
at r/a = 0.7. 



the other hand, only very recently it has been possible to investigate the frequency spectra 
also by mean of comprehensive nonlinear simulations. Indeed, long enough simulations are 
needed for such an analysis. Ref. |53| provides a first important example; nevertheless in this 
latter work, where the BES technique is used, the turbulence measurements are not resolved 
in kg. Here we report some preliminary results about the kg resolved frequency spectra and 
the comparison between measurements and nonlinear local GYRO simulations. 

The Doppler reflectometry technique measures the instantaneous spatial Fourier analysis 
of the density fluctuations, (5n (k, i) = Jy dx n {x,t) e^^'^. In the last paragraph we have 
discussed the information on the density fluctuations related to the wave-vector k, providing 
the kg fluctuation spectrum. On the other hand, thanks to the good time resolution of the 
diagnostic, a simple Fourier analysis on the same scattered signal can provide localized mea- 
surements of the kg resolved Sn frequency spectra. 

The Doppler shifted component corresponding to the scattering close to the cut-off layer, 
typically dominates these frequency spectra (or at least it can be clearly separated from any 
spurious signal that appears mainly at zero frequency). The parasitic components at / = 
can be due to the direct reflection of the beam, or of the refracted beam on inner elements of 
the chamber, or again to a backscattering all along the beam path. In most cases, however, 
the frequency spectrum clearly exhibits the Doppler backscattered component without any 
/ = spurious component, as detailed in Ref. 0S], revealing the good selectivity of the 
diagnostic. 

The interest in studying the fluctuation frequency spectrum for a given toroidal wave- 
number is strictly linked to the issue of the quasi-linear modeling. The argument has been 
already introduced in the paragraph |3.1.2| The statistical DIA theory for example, foresees 
a spectrum which presents a Lorentzian shaped broadening: this is also the choice adopted 
in our quasi-linear transport model QuaLiKiz. From the experimental point of view, the 
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Doppler measurements are generally not compatible with a such Lorentzian shape, neither 
with a Gaussian one, as detailed in Ref. [IS]. Therefore, a different approach has been 
developed for the interpretation of the experimental fluctuation frequency spectra, that will 
be here described. 

Since the 6n frequency spectrum corresponds to the fluctuation time-correlation function 
in the Fourier space, a classical statistical argument can be introduced as a very simple model. 
The approach closely recalls the discussion of the paragraph |3.1.1[ where a Lagrangian time 



correlation function (r) has been introduced in Eq. (3.15). Very simply, the mean square 
value of the displacement of a particle with a velocity v in the turbulent field, (A^), can be 
classically evaluated for times either long or short with respect to the Lagrangian velocity 
correlation time tl- This would lead to two limits, respectively (A^) = 2Dt = u^tlt and 
^A^) — v?t'^, where D is a diffusion coefficient and — {v^)- It clearly appears that the 
first expression corresponds in fact to a diffusive behavior, while the second to a convective 
one. The former considerations can be applied to a simplified model which describes a general 
transition from the convective to the diffusive behavior, leading to the expression: 

(A2) ("--l + e-^) (4.6) 



which correctly recovers the previous limits. 

An expression for the correlation function of the scattered field emitted by the Doppler 
reflect ometer's antenna has to be derived. This has been done in [321110]) giving the following 
expression: 

Ck(r) = /^e''"('''(*)-''^(*+"»\ ^5(k)F(k,T) (4.7) 
\ »j / t 

i.e. the product of the static form factor S (k) by the characteristic function of the probability 
of turbulent displacement F(k, r). The first one is a static factor which defines the kg 
spectrum previously considered, while the second one corresponds to the probability of the 
displacements. The latter one can be rewritten when assuming a given statistics for the 
distribution of these turbulent displacements (here denoted by A^-); under the (rather strong) 
hypothesis of normal statistics one has: 

i^(k,r) = (e''*"^-) P(A|T)e'''-^dA^ e''='(^'>/2 (4.8) 

Here the game is to find an explicit expression for the mean squared turbulent displacement 
(A'^). a general expression follows from the classical statistical approach of the Lagrangian 
dynamics of a particle with the velocity v |77j : 

{A^)^2u^f {T-s)Cy{s)ds with C„ (s) = (v (0) w (s)) /w^ (4.9) 
Jo 

where Cy is the Lagrangian velocity correlation function and = {v^^ is the variance of 
the velocity. A Lagrangian velocity correlation time is then given by tl = /q°° (t) dr. 
Under these hypotheses and using Eq. (4.6), the signal correlation function (^(k, t) results 
proportional to the following expression: 

C(k,r) oc (e^''-^^) = e-'=^(^^)/2 = e~'^^'<^~'^''^'^') (4.10) 
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In (4.10), the passage between the second and the third term has been possible thanks to 
the hypothesis of Gaussian statistics for the turbulent displacements; a rigorous justification 
can be found in |64j . making use of the mathematical theory of cumulants. In this work, Eq. 



(4.10) is here referred as the model test T function. 



In the Fourier space, Eq. (4.10) corresponds to a frequency spectrum at a given wave- 
number k, which can not be analytically calculated. Nevertheless, it is particularly useful to 
examine the frequency spectral shape resulting from the Fourier transform of the Eq. |4.10| 
in two relevant limits formerly cited: 

1. T ^ T^: Diffusive limit. The diffusive behavior leads to a Lorentzian broadening of 
the frequency spectrum. This was the expected result, since the similar argument of a 
diffusive random-walk of the particles in the turbulent field, produces the Lorentzian 
spectrum in the quasi-linear theory framework. Consequently, the wave-number k- 
dependence of the spectral widths FWHMk = 2Awk follows FWHMk oc k^. 

C'(fc,T) « e^*'^"^'^'^^ <^ 5fc (w) « Lorentzian ^ Aw^ oc fc^ 

2. T <C tl'. Convective limit. The convective component originates a Gaussian shape 
for the frequency spectrum. Differently from the Lorentzian broadening, this behavior 
results in a fc-dependence for the frequency widths of the type FWHMk oc k-'^ . 

C'(fc,T) sa e"''^"^'^^ 4^ S'fc (w) « Gaussian -i^ Aw^ oc fc^ 

Nevertheless, it has to be stressed that the previous expectations on the fc-dependence of 
the frequency broadenings are given in the hypothesis that both and tl do not depend 
on the wave-number k. Even if this could appear a strong limit, when assuming dominant 
and across the fc-spectrum, the previous scalings are meaningful to indicate dynamics 
dominated by diffusive or convective behaviors. 

The frequency spectrum corresponding to the whole T model, i.e. the Fourier transform 



of Eq. (4.10), presents a spectral shape which is intermediate between a Lorentzian and a 
Gaussian. Hence in this case, it is reasonable to expect that the /^-dependence of the frequency 
broadenings obeys to Aw oc k°'^''\ where the exponent has an intrinsic k dependence a — 
a (k), and 1 < a (fc) < 2. Under the same former hypothesis of and independent of k, 



the spectral exponent a can be numerically studied from the Fourier analysis of Eq. (4.10), 
recognizing a transition between the diffusive and the convective regimes. The properties of 
this simple T model can then be summarized as: 



C(fc,T)we V"^- J ^ S-fe (w) w intermediate Aw^ oc fc"^'^) 

The study of the fluctuation frequency spectra is here addressed comparing the Doppler 
reflectometry measurements to the nonlinear gyrokinetic simulations using GYRO. At this 
scope, the following relation is adopted to reproduce the experimental signals: 



Doppler 



(4.11) 



where T stands for a time Fourier transform. The latter expression is applied to the local 
nonlinear GYRO simulations done on the discharge TS39596. The time Fourier transform 
is performed over a sufficiently long interval of the nonlinear saturation regime predicted by 
the code. The density fluctuations are computed at the outboard midplane, according to the 
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diagnostic line of sight, as already described in the previous paragraphs. The wave- number 
k here refers to the poloidal wave-number ke selected by the Doppler reflectometry. Finally, 
since the GYRO simulations are done in the local limit, hence without any variation of the 
radial profiles, the average over the radial coordinate (. . .)^ is simply useful to cumulate 
statistics on the numerical predictions. This study refers to the same Tore Supra discharge 
already described, TS39596; also, the same radial position used for the k spectra analysis, 
r/a — 0.7, is considered. We refer to Tables 4.1 and 4.2 for the experimental parameters. 



The gyrokinetic nonlinear simulations are here used for the first time in order to quanti- 




Figure 4.10: (a) Frequency spectrum of Snj^ (where kgps — 0.82) from the GYRO simulation 
of TS39596 at r/a = 0.7 and comparison with different spectral shapes, (b) Experimental 
frequency spectrum from the Doppler refiectometry on the same discharge and radial position; 
best fits from Gaussian, Lorenzian and T model spectral shapes are also shown. 



tatively investigate the issue of the spectral shape of the frequency spectrum and to compare 
it to the experimental evidences. In Fig. |4.10[ the frequency spectrum of the density fluctu- 
ations at a given wave-number, computed by the GYRO simulation according to Eq. (4.11 ), 
is shown. The GYRO density fluctuation frequency spectrum is compared with a Lorentzian 
and a Gaussian spectral shapes. It clearly appears that both functions fail in reproducing 
the broadening observed in the nonlinear simulations. Nevertheless, it is worth noting that 
the Lorentzian and the Gaussian shapes can still be useful to obtain reasonable estimates on 
the values of the FWHM of the spectra. 

The T model described by Eq. (4.10) is moreover applied to the analysis of the GYRO 
results. Since there is not an analytical expression for the Fourier transform of Eq. (4.101, a 
nonlinear bisquare regression has been applied in the direct space, i.e. to the time correlation 
function of the density fluctuations predicted by GYRO; hence an inverse FFT is used to 
compare the frequency spectrum of the T model to the one directly obtain by the simulation. 
Surprisingly, the T model shows an extremely good agreement in reproducing the GYRO 
frequency spectral shape (a factor = 0.9996 is obtained from the nonlinear regression). 
This result is even more relevant when compared to the experimental evidences. Until now, 
the T model is often used as practical solution for the fitting of the experimental frequency 
spectra of the Doppler reflectometry [48 , but without a clear understanding of the reasons 
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of its superiority with respect to the Lorentzian and the Gaussian shapes. In Fig. 4.10[ b), 
the experimental Sn frequency spectrum at a given wave-number kg selected by the Doppler 
reflectometry, is shown for the same discharge TS39596 at r/a = 0.7. The asymmetry seen 
on the positive side of the spectrum is most probably due to spurious components coming 
from the backscattered signal all along the beam path. The experimental spectral shape ex- 
hibits analogous properties with respect to the GYRO results, with a frequency broadening 
intermediate between a Gaussian and a Lorentzian function, but conversely well fitted by 
the simple T model. The present study provides then a very good validation of the methods 
used in the analysis of the Doppler reflectometry. 

Following from these results, the consequent crucial test for the validation of the nonlin- 




Figure 4.11: kg dependence of the FWHM relative to the Sn frequency spectra on TS39596 at 
r/a = 0.7: the measurements by the Doppler reflectometry are compared to the expectations 
by GYRO and from the T model. 



ear simulations, is finally the quantitative comparison of the fc-dependence of the frequency 
broadenings. Previous experimental studies addressed this issue using a GO2 laser back- 
scattering technique for the analysis of the density fluctuations. The findings revealed that 
the frequency widths Afk increase at higher k according to Afk (x k", where 1.1 < a < 1.5. 
More precisely, a transition of the spectral exponent was observed across fc: a ~ 1.5 was 
found at low fc, while decreasing for large k, a ss 1.1. 

On the discharge TS39596, the radially localized Doppler measurements at r/a = 0.7 allow 
reliable evaluation of the FWHM of the Sn frequency spectra associated to each wave-number 
kg. These are quantitatively compared to the FHWM computed by the GYRO simulation: 
the results are shown in Fig. |4.11| 

Firstly, the GYRO simulation shows a remarkably good quantitative agreement with the 
experimental spectral widths from the Doppler reflectometry. The experimental data are 
characterized by a non negligible dispersion, but they can be fitted by a power law whose 
mean spectral exponent is aooppier = 1.4 ± 0.4. The predictions coming from the simple 
T model are in remarkable agreement with both the measurements and the nonlinear gy- 
rokinetic simulation. More in particular, the T model fc"^'') dependence for the broadenings 
exhibits a slight decrease of the spectral exponent a (k) at larger fc, as seen by the nonlinear 
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simulation. These evidences provide a relevant additional information with respect to the 
former result of Fig. |4.10[ obtained at a given wave-number kg. Here in fact, the whole 
fc-dependence of the fluctuation frequency spectra is consistent with the outcomes of the T 
model. 

Summary 

In summary, the turbulent dynamics studied by mean of both density fluctuation mea- 
surements with the Doppler reflectometry and nonlinear gyrokinetic simulations, appears in 
good quantitative agreement with the following main hypotheses of the simple T model: 

1. In the nonlinear saturation phase, the particle displacements in the turbulent field (A^) 
are well described by a simple transition between diffusive and convective behaviors. 



described by Eq. (4.6). Individually, neither the first nor the second mechanism are 



able to reproduce the observed frequency spectral shape. 

2. The distribution of the particle turbulent displacements appears to be consistent with 
the hypothesis of a normal Gaussian statistics. 

3. Assuming dominant and tl as fitting parameters independent of the wave-number 
k, the A:"''^' dependence of the frequency broadenings is correctly reproduced for a wide 
spectrum of wave- numbers k. 

Despite the success of this very simple model in reproducing both the Doppler measurements 
and the nonlinear gyrokinetic simulations, the issue of the fluctuation frequency spectra 
certainly deserves additional study in the future. In particular, the analysis with more so- 
phisticated theoretical models, namely the advanced statistical turbulence theories, is highly 
desirable for relevant advances in this research area. 



4.2 Modeling of the nonlinear saturation process 

In the last section we have treated the validation of the nonlinear gyrokinetic simulations 
against turbulence measurements, concerning both scalar, i.e. rms Sn/n (r) and Xef f and 
spectral quantities, i.e. Sn/n (fcg, k^) and Sn/n (fc,a;). The aim of the present section is then 
using the information carried by these nonlinear simulations in order to produce accurate 
hypotheses for the saturation model required by the quasi-linear modeling. In order words 
the question is: how to model the saturation spectral intensity 

It is important to clarify that the answer to the latter question can not be immediately 
derived from the direct comparison between the turbulence measurements and the nonlinear 
simulations illustrated in the previous section. The reason is simply due to the fact that 
the measured quantities do not coincide with what is required by the quasi-linear transport 



modeling. In fact, recalling the quasi-linear expressions Eqs. (3.19 1-(3.20 1, it appears that 
(1) the spectral saturation intensity in the quasi-linear formulation refers to the potential 
fluctuations, and not to the density fluctuations measured by the diagnostic^ (2) \S(j)k\'^ has 
to be calculated according to a proper flux-surface average, and not at the given position in 
the poloidal cross section where the diagnostics measure the fluctuations (typically on the 
low field side). 



^The approximation Sn/n e5(f>/Te is in fact rigoursly valid only if neglecting the non-adiabatic electron 
response 
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The model of the saturation intensity ^l^^/cl^y used in QuaLiKiz can be written according 
to the expression: 

'\5<j)u?)=SkY,Lk,j (4.12) 

3 



with the explicit summation over the j hnear unstable modes. The expression (4.12) is 
composed of two different parts: 

1. The first term Sk, is the spectral functional shape in the k space, which is assumed the 
same for each unstable mode. 

2. The second part ifcj-, concerns instead the saturation level that weights the j unstable 
solution at each wave-number k. 

After a brief discussion on the issue of the frequency spectrum in the actual version of Qua- 
LiKiz, the following two paragraphs will be dedicated respectively to these two issues. 



4.2.1 Choices on the frequency spectrum 

At this point it is important to mention that the issue of the frequency spectrum, just dis- 



cussed in the previous paragraph, does not explicitly appear in the expression (4.121. In 



fact in QuaLiKiz, as shown in the paragraph 3.1.2, the choice of a particular frequency spec- 



trum for the saturated potential is included in the formulation of the quasi-linear response. 



through the additional imaginary contribution in the resonant denominator (Eqs. (3.271 



(3.28)). Namely, in QuaLiKiz, the frequency spectrum is assumed with a Lorentzian shape 
characterized by broadening equal to linear growth rate, i.e. a choice which is not consis- 
tent with the evidences of the nonlinear gyrokinetic simulations, as shown by Figs. 4.10||4.9 



Therefore this point clearly deserves dedicated further studies and improvements with re- 
spect to the present formulation in QuaLiKiz have to be done. Nevertheless interestingly, as 
it will be shown in the next chapter, the actual choices on the frequency spectrum provide a 
good agreement on the total turbulent fluxes by QuaLiKiz with the expectations by nonlinear 
gyrokinetic simulations. 



4.2.2 Spectral structure of the saturated potential 

Local nonlinear gyrokinetic simulations are again used in order to determine the k spectral 
shape of the saturated potential required by the quasi-linear model. It is important to remind 
that nonlinearly, the k spectrum of the saturated potential can not be separated into distinct 
contributions corresponding to the linear modes j, as already pointed out in the paragraph 
|3.2.3| For this reason in QuaLiKiz, the k saturation spectral shape is assumed identical for 
all the linear unstable modes j. This point remains an arbitrary hypothesis that could be 
object of further studies. 

The quantity that has to be computed from the simulations is (\^5(f>{kg)\^^ , where (. . .) 

corresponds here to the time and flux-surface averages. The result from the GYRO simulation 



of TS39596 at r/a = 0.7 are reported in Fig. 4.12 



Even if for the purposes of the quasi-linear modeling only the dependence of the saturated 
potential on the toroidal wave-number n, hence the kg spectrum, is important, also the kj. 



spectrum is shown in Fig. 4.12 for comparison. It appears that the power law k^p^ fits 
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Figure 4.12: kg and kj. spectra of the fluctuating potential from the GYRO simulation of 
TS39596 at r/a ^ 0.7. 



very well the nonlinear results in the range 0.2 < kePs < 1-0, i.e. for the kg scales above the 



peak corresponding to fc"' 



6. max ' 



Interestingly, the kpj^ cascade found by the nonlinear gyrokinetic simulation of Fig. 4.12 
is well reproduced also by a recently proposed analytical model for the k spectrum [44; . This 
work deals with the minimal hypotheses of homogeneous and isotropic fluctuations, dealing 
with a weak turbulence theory, in particular the nearly adiabatic limit of the Hasegawa- 
Wakatani model [46]. In the limit of dominant nonlocal disparate scales interaction through 
the large scale flows (zonal flows), the following analytical expectation is derived: 

One of the most original feature of this model is the presence of a change of slope in the 
nonlinear k cascade around kps « 1, leading to a steeper spectral exponent for the smaller 
spatial scales kps > 1.0. Even if this point has been found in nice agreement with laser 
back-scattering turbulence measurements [37], it will not be examined here, since the higher 
k components of the spectrum kps > 1.0, are not resolved by our nonlinear simulations, 
because of the computational cost. On the other hand, in the range 0.2 < kps < 1.0, both 
the analytical model and the nonlinear simulations coherently find the spectral shape k±pj^. 

In the light of these results, in QuaLiKiz, the k saturation spectral shape Sk is chosen 
according to the following relations: 



Sk oc kpl kps < kf„^^^ps 
Sk oc kp-"^ kps > k'^]„^^^p, (4.14) 



The symmetry of the spectrum model around kg'-^^^ps described by (4.14), is justified by 



both experimental results using the beam emission spectroscopy [74] and nonlinear simula- 
tions [52]. 



4.2 Modeling of the nonlinear saturation process 
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4.2.3 Saturation rules 



In order to corapletely define the saturated k spectral shape Sk adopted in the quasi-linear 
calculation, a choice for the value of the nonlinear peak k'^^^^x appearing in Eq. (4.141 has 
to be made. Generally, the nonlinear spectrum peaks in the kg axis for values keps ~ 0.2 
(see Fig. 4.12 1, i.e. lower than the typical linear growth rate maximum of the most unstable 

and 



mode kgps ~ 0.4. This feature has been confirmed by both experimental evidences ^ 

numerical simulations [17]. Moreover, a relevant dependence of the value of fc^'^^^ on the 
safety factor q is observed in nonlinear simulations [271 151) . In QuaLiKiz, a simple mixing 
length argument, as introduced in the paragraph 3.1.1 Eq. (3.18), is used with two different 



goals: firstly, to produce a pertinent choice of kg^max' and secondly to provide the saturation 
level weighting the contributions from the different linear unstable modes. 

The value of kg^^^^ is chosen such that the mixing length diffusivity factor is maximum, 
according to: 



max(I?e//) 



7 



(4.15) 



In Eq. (4.15), both the growth rate and the (k^) mode structure refers to the linear most 
unstable mode, so that the resulting kg^^^^ unambiguously defines a saturation shape Sk 
which is kept identical for all the linear unstable modes. 

In QuaLiKiz, (k^) is chosen according to the form originally proposed by [55l[27l|4], valid 
for strongly ballooned modes and adding the impact of the MHD parameter a: 



where 



in 



J 9^ |(5(/)5^™ {e)\ de 
{e)fde 



I\ 



(4.16) 



(4.17) 



In Eq. (4.17), the 6 ballooning structure of the linear eigenmodes is involved; in QuaLiKiz 
a trial Gaussian eigenfuction is used, such that (50'*" [0) oc y/we~ ™ , where w is the 

mode width, solution of the fluid limit, and d is the distance between two resonant surfaces. 
In that case, Eq. (4.171 then reduces to: 



2d^r(3/4) 

w2r(i/4) 



(4.18) 



where T is the mathematical Gamma function. 

The functional shape of the saturation k spectrum Sk of Eq. (4.14) can be specified 
with the former hypotheses derived from a mixing length argument. Finally, the last choice 
refers then to the absolute value of Sk- It is reasonable to expect that the linear unstable 
eigenmodes contribute in a different way to the total turbulent flux, according to the dominant 
or sub-dominant nature of the mode. Hence, the saturation model adopted in the quasi-linear 
transport modeling has then to account for weighting factors on the different linear unstable 
roots. In the actual version of QuaLiKiz, each unstable mode j is weighted through the 
mixing length hypothesis using its corresponding growth rate and mode structure. In other 
words, recalling Eq. (4.12), we have: 
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This choice naturally implies that the largest contribution to the total turbulent flux is car- 
ried by the leading linear mode at each wave-number kg, while the sub-dominant solutions 
bring smaller terms. Nevertheless, these sub-dominant terms can have a crucial role in con- 
ditions where two modes are competing with comparable growth rates, while retaining only 
the most unstable solution would imply a failure with respect to the nonlinear expectation 
for the total turbulent flux. 

The weighting factor adopted in QuaLiKiz expressed by Eq. (4.19) is completely equivalent 
to the choice adopted in the TGLF model, demonstrating that this rule is able to achieve a 
good agreement with the predictions by nonlinear gyrokinetic simulations. Analogous good 
agreement will be treated in the next chapter. 



Summary 



Finally, the fundamental choices on the saturation model adopted in QuaLiKiz are here 
summarized: 

1. A general functional shape for the k saturation spectrum Sk is identified, characterized 
by a nonlinear cascade with a spectral exponent a = —3 and symmetric around the 
nonlinear peak kg^^^^ . The validity of this model is coherently supported by nonlinear 
gyrokinetic simulations, turbulence measurements and simplified analytical models. 

2. The choice of a Lorentzian frequency spectrum with a finite broadening equal to the 
linear growth rate at each wave-number k is assumed in the formulation of the quasi- 
linear response. This hypothesis does not agree with the results inferred by both 
nonlinear simulations and turbulence measurements. Actually, further improvements 
on this point can be investigated in future work. 

3. A mixing length argument for the effective turbulent diffusivity is used to derive both 
the nonlinear peaking kg^^^^ in the model of the k spectrum, and the weighting factors 
for the different linear unstable modes. 



Chapter 5 

Operating the quasi-linear 
model 

In the previous chapters, both the vahdity of the gyrokinetic quasi-hnear response and the 
improvement of the model of the saturated potential have been accurately discussed. Finally, 
in this last chapter, the two main parts of the quasi-linear transport model are put together, 
in order to produce the estimates of the total turbulent fluxes of energy and particles. Two 
levels of validation will be here discussed. 

The first one is the direct comparison between the quasi-linear turbulent fluxes and the 
predictions by the nonlinear gyrokinetic simulations. The first section will be then devoted to 
the verification of the quasi-linear QuaLiKiz versus the nonlinear GYRO expectations. Sev- 
eral parametric scans are presented for tokamak relevant conditions, allowing to investigate 
in which conditions the quasi-linear approximation is able to track the nonlinear results. 

The second level of validation concerns the application of the quasi-linear transport model 
to realistic tokamak scenarios. The aim is here gaining interpretative and predictive ca- 
pabilities on the actual experimental observations. This step necessarily implies that the 
quasi-linear model is coupled within an integrated transport solver. The latter integrates the 
quasi-linear turbulent fluxes with the sources of particles and energy and the reconstruction 
of the magnetic equilibrium, solving the time evolution of the plasma profiles. 

The second section will be dedicated to the coupling between QuaLiKiz and the inte- 
grated transport code CRONOS, while the third section deals with some first analysis of the 
experimental results. 



5.1 Validation against the nonlinear predictions 
5.1.1 Parametric scans 

In this paragraph the quasi-linear turbulent fluxes computed by the actual version of Qua- 
LiKiz are compared to local nonlinear gyrokinetic simulations using GYRO. GYRO has been 
run using the standard numerical resolution summarized in the paragraph |3. 2. 1[ The funda- 
mental quantities that are compared are the effective energy and particle diffusivities, defined 
according to the relations Fj — —Ds^rng and Qs — —UsXs^rTs for each species s. Additional 
details on the calculation of the total turbulent flux by the nonlinear GYRO simulations are 
given in Appendix [B| 



The quasi-linear fluxes computed by QuaLiKiz follows from the Eqs. (3.19 )-(3.20 1, using 
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the model for the saturated potential described in the previous chapter |4] Finally, a single 
arbitrary constant Co is used to renormalize the quasi-linear results: in other words, Co is the 
factor responsible of rescaling the quasi-linear over-prediction with respect to the nonlinear 



results, as discussed in the paragraph 3.2.3 Practically for QuaLiKiz, Co = 1.6~^ has been 
chosen in order to renormalize the value of the quasi-linear ion energy flux to the GYRO 
nonlinear prediction for the GA-standard case. This value is fixed and applies to all the trans- 
port channels (energy and particle) and to all the plasma parameters used in the simulations. 




Figure 5.1: Ion energy, electron energy and particle effective diffusivities from GYRO (points) 
and QuaLiKiz (lines) for the R/Lt scan based on the GA-standard case. 



In the first scan, both the ion and electron temperature gradients are simultaneously 
varied on a wide range: 4.5 < R/Lti = R/Ltc < 13.5. The effective energy and particle 
diffusivities are expressed in GyroBohm units xgb = P^Cg/a. Across the whole scan, the 
ion and electron energy and particle fluxes computed by QuaLiKiz and GYRO agree within 
15%. Both the ratio between the transport channels and the parametric dependence are well 
captured by the quasi-linear approach. 



Focusing on the R/Lt = 9.0 point of the scan presented in Fig. 5.1 i.e. the GA-standard 



case, a first preliminary test is done exploring the coupled ITG-TEM-ETG transport. Nowa- 
days, this kind of exercise is still challenging beacause of the extremely high computational 
resources needed to run nonlinear gyrokinetic simulations retaining disparate spatial scales. 
Here we refer to the results by the first coupled ITG-TEM-ETG gyrokinetic nonlinear sim- 
ulations reported in |100| . using the GYRO code. In this work, both the medium-large ion 
and the small electron spatial scales up to kgpe 1 are coherently solved by the nonlinear 
simulations, but using a reduced mass ratio \/mi/ me = 20 instead of \fm^rne — 60. 
From the point of view of the quasi-linear modeling with QuaLiKiz, the actual hypotheses 
described in the previous chapters are kept without any modification. In particular it is 
worth to stress that the spectral exponent of the fc saturation spectrum is kept a = —3, 
while there are experimental and theoretical insights suggesting a steepening of the spectral 
slope for sub-ion spatial scales kgpi > 1, as discussed in the paragraph |4.2.2 Referring to 
the GA-standard case, QuaLiKiz predicts that the ETG contribution (i.e. for kgps > 1) to 
the total electron energy flux is 11%, in very good agreement with the value of 10% obtained 
by the massive GYRO simulation. 



5.1 Validation against the nonlinear predictions 
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As important point it has to be noted that the GA-standard case does not contain any ad- 
ditional E X B shearing rate coming from an external radial electric field. The latter one 
has been found in Ref. [lOOj . to be very efficient in quenching the low k ITG-TEM trans- 
port, while leaving almost unaffected the ETG contribution. This evidence has suggested 
[59] that in certain experimental conditions, with high E x B shearing rates, the ETG driven 
transport can still become dominant, since the ITG-TEM contributions are strongly reduced. 

A second example is a direct application to an experimental collisionality Vi, scan realized 
on Tore Supra plasmas. This has been realized as a dimensionless scaling experiment, since 
both the other two relevant dimensionless parameters, ps and /3 are kept constant across 
the different discharges, thanks to a coherent variation of the magnetic field and the tem- 
peratures, so that only the collisionality is changed. The main plasma parameters are here 
summarized in Table l5Tl 

The scaling of transport is a crucial test for quasi-linear models. Two effects are 



Ro/a 


r/a 


R/Lxi 


R/Ltc 


R/Ln 


q 


s 






p* 


/3 


3.25 


0.5 


8.0 


6.5 


2.5 


1.48 


0.72 


1.0 


1.0 


0.002 






Table 5.1: Plasma parameters of the Tore Supra dimensionless collisionality scan at r/a — 
0.5. With respect to the experimental values, only /? is artificially set to in the GYRO 
simulations. 




Figure 5.2: Ion energy, electron energy and particle effective diffusivities from GYRO (points) 
and QuaLiKiz (lines) for the collisionality scan based on the Tore Supra discharges. 



potentially at play with opposite consequences on the total flux levels. The collisional damp- 
ing of sheared flows would in fact result into an increase of the turbulence level at higher 
collisionality; hence the fluxes would be enhanced, as pointed in Ref. [3S]. This effect is not 
taken into account in QuaLiKiz, which does not include sheared flows. Conversely, the linear 
collisional TEM damping would act in the opposite direction, reducing the linear instability 
drive as therefore the total turbulent flux at higher collisionality. 

For the plasma parameters here explored, the GYRO simulations find a visible reduction of 
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the transport level when increasing the coUisionality, as shown in Fig. |5.2[ Therefore the 
linear TEM damping is found to be dominant in these nonlinear simulations with respect to 
the coUisional drag of sheared flows. Finally, Fig. |5.2| demonstrates that, for experimental 
values of coUisionality, the quasi-linear modeling by QuaLiKiz is able to well reproduce the 
nonlinear diffusivities predicted by GYRO simulations, performed with pitch-angle scatter- 
ing operators on both electrons and ions. The coupled dynamics between ion and electron 
non-adiabatic responses is crucial for both GYRO and QuaLiKiz. In particular, the particle 
flux reverses direction as ;/* increases, as already noticed in Ref. [4J. 

The third scan illustrates a variation of the ratio Ti/T^ performed on the GA-standard 
case. Both the nonlinear and the quasi-linear results are reported in Fig. |5.3| 
Some discrepancy between the quasi-linear particle flux predicted by QuaLiKiz and the 




Figure 5.3: Ion energy, electron energy and particle effective diffusivities from GYRO (points) 
and QuaLiKiz (lines) for the T^/Te scan based on the GA-standard case. 



nonlinear result by GYRO is found for Tj/Tg — 0.5. For this point of the scan in fact, ITG 
are expected to be dominant, because of the proportionality of the linear ITG threshold with 
Ti/Tg. The disagreement between the quasi-linear and the nonlinear particle flux is then 
coherent with the previous results, that indicated the weakness of the quasi-linear modeling 
when dealing with inward ITG particle flows (see paragraph 3.2.3). Actually this aspect 
appears as one of the most significant failures of quasi-linear modeling. 



Fig. |5.4| illustrates the TEM to ITG transition on the GA-standard case, realized keeping 
fixed R/ Ltc = 9.0 and varying only the ion gradient R/ Lti- 

The electron energy fluxes are well matched; discrepancies are instead observed on the 
particle fluxes for strong ITG turbulence and for the ion energy flux for R/Lxi < 4 and 
R/Lxi > 13. At R/Lxi < 4, TEM become the dominant unstable modes, while the marginal 
conditions for ITG turbulence could be responsible for this quasi-linear failure on the ion 
energy flux. The quasi-linear model TGLF for example, adopts a modifled mixing length 
rule that is proportional to 7 + 7^ '^ instead of the simple linear relation with 7. On the 
other hand, above R/Lti = 13, the QuaLiKiz overestimations can be ascribed to a more pro- 
nounced effect of zonal flows in the nonlinear saturation of the ion energy transport for ITG 
dominated turbulence. Hence, the impact of both sheared flows and marginal turbulence on 
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Figure 5.4: Ion energy, electron energy and particle effective diffusivities from GYRO (points) 
and QuaLiKiz (lines) for the TEM to ITG transition based on the GA-standard case, fixing 
R/Ltc = 9.0 and varying R/Lxi- 



the quasi-linear transport modeling deserves dedicated additional analysis. 

Fig. |5.5| presents a dilution scan operated on plasmas with D main ions, electrons and He 
impurity on the GA-standard case. Moreover it is assumed that: Tf/e = ^i, R/ Lthc — R/Lti 
and R/LnHe = R/Ln- For both energy and particle transport of all three species, the 
quasi-linear predictions by QuaLiKiz are shown to rather well agree with GYRO simulations 
realized with full nonlinear gyrokinetic treatment of the three plasma species. 




Figure 5.5: Effective energy and particle diffusivities from GYRO (points) and QuaLiKiz 
(lines) for a dilution scan, based the on GA-standard case considering D ions and electrons 
plasma with He impurity. 
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5.2 Coupling with an integrated transport solver 

This section is dedicated to the coupling of the gyrokinetic quasi-linear transport model 
QuaLiKiz into a system of time dependent fluid transport equations for the evolution of 
the plasma profiles. Starting from a general formulation, the discussion will be particularly 
referred to the application within the integrated transport code CRONOS. 



5.2.1 The time dependent transport equations 

The definition of the fundamental quantities which are evolved in the transport equation 
directly follows from the velocity space integration of the kinetic Fokker-Planck equation 
for the species s, which describes the evolution of the distribution function fs (x,v,t). The 
velocity moments from the 0*^ to the 3'''^ order are respectively the density, the particle flow, 
the stress tensor (pressure) and the energy flux: 

n^fd^vf T = nu=[(fvvf (5.1) 



P = pi + 11 = m / d^v V X v/ 



Q = I d^v -mv'^vf 



(5.2) 



where the scalar pressure p is defined from the stress tensor as Tr (^P^ /3. Another important 

3'^'^ velocity moment is the heat flux, i.e. the energy fiux in the frame of the moving fluid. 
Two alternative definitions can be found in literature 1501 : 



q = y d^v -m\v — u]"^ (v — u) f <l' ~ J d^ v -mv'^ (v — u) f 



(5.3) 



Using the hypothesis that the distribution function / can be decomposed into a Maxwellian 
component plus an expansion based on the parameter e — p/L, i.e. the ratio between the ion 
gyroradius and the typical macroscopic transport length, at the 0*^ order in e Eqs. (5.3 1 are 
linked to the other fluid moments according to: 



Q 



Q 



(5.4) 



At this point it is particularly relevant to make a link with the quantities that are com- 
puted by QuaLiKiz. The latter one in fact deals only with the fluctuations of the distri- 
bution function Sf with respect to a local Maxwellian equilibrium. Recalling the QuaLiKiz 
expressions (3.19)-(3.20) for the particle and energy turbulent flows, these corresponds to the 



quasi-linear estimates of the following quantities: 



r«' w J d^v w6f ■ Vp ~ J 



d^v -mv^vSf ■ Vp 



(5.5) 



Eqs. (5.51 have to be 
(5.1|-(5.2). It appears 



where p represents the perpendicular cross-fleld radial coordinate, 
compared to cross-fleld components of the flows expressed by Eqs. 
that the these two approaches are compatible only if assuming that the neoclassical and the 
turbulent transport are additive. The coUisional neoclassical transport can in fact not be 
captured by the 6f approach used in the quasi-linear modeling. In other terms, a possible 
influence of the turbulence on the evaluation of the neoclassical terms is here neglected. This 
point represents still nowadays an open subject of research and the validity of the hypothesis 
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of additivity remains to be addressed. 

The fundamental time dependent transport equations stem from the conservation of parti- 
cles and energy in the system. These conservation equations can be derived taken respectively 
the 0*** and the 2"^^ velocity moments of the Fokker-Planck equation. The conservation of 
particles is then expressed by: 

^ + V • (nu) - S (5.6) 
ot 

where the RHS represents a source term, while the equation for the energy conservation 
reads: 

^^+V-Q = Qc + u-(F + enE) (5.7) 



In Eq. (5.7), the two terms Qc and F result from the integration of the collision operator 
C (/), defining the coUisional energy exchange Qc = J (Pv ^ |v — u| C (/) and the friction 
force F = m J (Pv v x v C (/). Collisions have to respect the conservation of momentum and 
energy, so that it can be written: 

^F, = ^Qc. +u, -F, =0 (5.8) 

s s 

The main goal of the integrated transport modeling codes is to solve the flux-surface av- 



eraged version of the local continuity equations (5.6)-(5.7l for the particles and the heat flux. 
The key point of this kind of approach is the reduction of the dimensionality of the problem: 
the original 6D plasma dynamics (3 dimensions in the real space and 3 dimensions in the 
velocity space) is reduced to a ID approach. The flrst reduction, lowering by a factor of 3 
the dimensionality, corresponds to the cPv velocity integration of the Fokker-Planck equation 
(fluid approximation). A further reduction of 2 dimensions is linked to the fact that, at the 
first order in the expansion parameter e, the thermodynamical plasma quantities such as 
density, temperature and pressure are constant over the flux surfaces of an axis-symmetric 
magnetic equilibrium. The latter one is usually numerically solved in a 2D space, but the 
transport processes can be simplified to a ID approach, which deals only with the radial 
cross-field transport. This procedure requires then a proper flux-surface averag^ 

The average of a scalar quantity z over a given magnetic surface can be expressed con- 
sidering an infinitesimal volume around the flux surface: 

In the case of the continuity equation for the electrons, this can be obtained from the flux- 



surface average of Eq. (5.6). The cross-field electron flow can be written as Fe = Ue (u^ — Up)- 
Vp, where Up represents the possible velocity of the flux-surface labeled by the coordinate p 
in the laboratory frame. The continuity equation equation then reads: 

|(^^W) + ^(r(r,)) = y'(5,) (5.10) 

Actually, the equation which is effectively solved by most of the integrated transport codes, 
including CRONOS, prefers to use as unknown, a modifled form the particle flux, expressed 



^The system solved by the integrated transport codes can also be referred as a 1.5D problem, since an 
axis-symmetric magnetic equilibrium is solved in a 2D space, but the transport processes are mapped into a 
ID grid. 
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The issue of the magnetic equihbrium geometry deserves some additional care. Firstly, 



the transport equations solved by CRONOS, like Eq. (5.101, employ flux-surface averaged 
quantities, while the quasi-linear turbulent fluxes computed by QuaLiKiz do not depend on 
the poloidal angle. Therefore, the metric factor (|Vp|) has to be introduced in order to 
account for the quantities averaged across a surface of constant p. This affects both the 
profile gradients provided as input by CRONOS to QuaLiKiz and the resulting quasi-linear 
estimates for the turbulent flows: referring for example to the case of the particle flux, one 
has: 

(Vpn) =apn(|Vp|) and (r*') = T*' (|Vp|) (5.11) 

Secondly, as here above explained, the CRONOS code has to possibility to treat a generally 
shaped axis-symmetric magnetic equilibrium, while QuaLiKiz solves the linear gyrokinetic 
dispersion relation on a simplified s — a circular equilibrium. Hence, the direct impact of the 
shape of the magnetic geometry on the characteristics of the turbulence is not captured by our 
quasi-linear model. Recently, this aspect has been object of several studies using nonlinear 
gyrokinetic simulations 1501 13] . Both global PIC and local Eulerian simulations recognized 
a stabilizing effect of the elongation on the ITG-TEM turbulence, in terms of a reduction 
of both the linear growth rates and the nonlinear saturation level. Even if these effects are 
presently completely neglected in QuaLiKiz, several strategies can be explored in the future 
in order to operate an heuristic correction. The TGLF model for example, uses an adjusted 
factor for the total quasi-linear flows, which allows to achieve the best flt with a large database 
of local GYRO simulations run with a shaped Miller magnetic equilibrium. Otherwise, in 
Ref. [3], it has been proposed an effective normalized gradient length (obtained by dividing 
the real gradients by the square root of the elongation) , which is shown to correctly reproduce 
the linear growth rates in presence of finite elongation. 



The solution of the continuity equation Eq. (5.10) requires some cautions. Pinch terms 
are very important for the the particle transport, since they usually balance the diffusive 
terms to produce a net flux close to zero. For this reason in CRONOS, the particle flux 



appearing in Eq. (5.10) is explicitly split into two separate diffusive and convective terms, 
according to: 

(Fe) - F: (|Vp|') = -D, (|Vp|') dpn, + K (|Vp|') He (5.12) 

The latter expression has to be compared to the quasi-linear flow calculated by QuaLiKiz, 
that can be computed separating the diffusive and convective components, and it can be 
written as: 



{Tf) 



-Df (V,n,) + {y,t + Vf) (5.13) 



(|Vp|) 

where Vth and are respectively the thermo-diffusion and the compressibility components 



of the convective velocity. Finally, recalling Ec^s. (5.111, the QuaLiKiz coefficients that are 
used in CRONOS in order to solve the continuity equation ( |5.10 1 are]^ 



D, = Df and V, = (v.t + vA ^^^^^ (5.14) 



^For the diffusion coofficiont , tfic approximation (|Vp| ) ~ is used. 
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The flux-surface averaged version of the electron heat transport equation follows from 
the general expression for the energy conservation Eq. (5.7). Recalling that CRONOS uses 
the conventional definition for the heat flux, i.e. q in Eqs. (5.3)-(5.4|, after several algebra 



(see |54j for details) one obtains the equation for the radial cross-field electron heat transport 
ge = Qe • Vp: 



2 V'i dt 

-V (u 



{iPe)V'i) 



d_ 

dp 



v 



2 ' ' 



Vpe) - V {Qc,^) - V (U,; • (F,; + Z,71,E)) + V {J ■ E) + V (Se) (5.15) 



It is worth noting that both the latter heat transport equation (5.15), and the previous 
continuity equation Eq. (5.10), are derived retaining a finite velocity of the flux-surface in 
the laboratory frame Up, while QuaLiKiz computes the turbulent fluxes only relative to the 
magnetic surface , i.e. to the velocity u — Up. That implies that the total electron heat flux 
defined as Qe = Qe — fpUg, following from Eq. (5.4 1, has to be approximated because of 
the missing Up term in the QuaLiKiz prediction. More precisely, the flux-surface averaged 



quantity appearing in the second term of Eq. (5.151 can be re- written as 



-Pe (Ue - Up) 



Vp ) = ( Q, 



Vp 



(5.16) 



The approximation can then be easily expressed when neglecting the (|peUp • Vp) term, 
i.e. this term reduces to the total energy flux (qe -I- ^TeTe) (Qe • Vp) Even if this point 
represents a limit of the actual formulation, in most practical cases, the corrections due to 
the movement of the magnetic surfaces are negligible. 

Likewise for the continuity equation, CRONOS solves Eq. (5.16) in terms of the variable 

9e 



(|Vp|^) 



On the other hand, the flux-averaged energy flux computed by QuaLiKiz reads 



as {Qf) = Qf (|Vp|). Using the just mentioned approximation, the QuaLiKiz term which 
enters in the CRONOS electron heat transport equation is the following 



(|vp| ) 

ivpr 



(5.17) 



A brief description of the RHS terms of the Eq. (5.151 is here given. The first one repre- 



sents an apparent energy flux due to a finite velocity Up of the flux-surfaces. As previously 
observed, this additional component is usually negligible in comparison to the other source 
terms, and in any case it is not retained by QuaLiKiz. The second one refers to the non- 
negligible coUisional energy exchange from the electrons to the ions. This terms is calculated 
by CRONOS according to the derivation of Ref. [SD]. The third term formally represents 
the friction forces between the electron and the ion species: this small correction is treated 
in CRONOS using the neoclassical theory. The last two terms correspond to the heating 
sources, respectively to the J • E ohmic heating induced by the plasma current, and to exter- 
nal heating sources, like the RF power and the Neutral Beam Injection. 



■^Practically, CRONOS solves the Eq. l |5.15| splitting the diffusive and convective terms, in order to avoid 
injecting a negative diffusion coefficient, analogously to the continuity equation for the particles. On the 
other hand, such explicit splitting is not done for the energy transport inside Qualikiz, in order to speed up 
the calculations and because the energy flux is almost always outward, i.e. it corresponds to Xeff > 0- For 
this reason the QuaLiKiz energy flux is initially stored as a purely diffusive effective transport coefficient, 
while a purely numerical (_D, V) splitting is done in the CRONOS solver in the case of of inward heat ffux. 
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The equation that governs the cross-field ion heat transport qi = 
analogous to the electron one. It reads: 



Vp is completely 



-V 



{l<Pr)V'i 



d_ 
dp 



(Up • Vpi) + V (ge,,;) + V (u, • (F, + Z,n,E)) + V {S,\ 



(5.18) 



The whole discussion following from the previous electron heat equation can be readily applied 
here. It is worth noting that the second and the third term of the RHS are identical to the 



Eq. (5.15), but with opposite sign. 



CRONOS only solves the Eq. (5.18) for the total ion pressure, assuming the same temperature 
for all the s ion species, i.e. Pi^tot = Ti'^^ns- Conversely, QuaLiKiz can presently account 
for electrons, ions and an additional ion species, computing separate particle and energy 
fluxes. With analogous notation with the previous one, the QuaLiKiz estimations that enter 
in the CRONOS total ion heat transport equation are finally: 



:T,r* 



irnp 



(|Vp| ) 
ivpr 



(5.19) 



5.3 Application to the experiments 
5.3.1 Preliminary results 

The present version of QuaLiKiz is being coupled and tested within the CRONOS code. 
Unfortunately, only preliminary results are here presented, since at the time of the redaction 
of the thesis this is still an ongoing work, which certainly deserves much more attention. An 
important task that is currently being realized is due to the fact that for each iteration of 
the transport solver, QuaLiKiz has to be run at each position of the radial grid. For this rea- 
son, in order to perform time-dependent simulations within a reasonable time, the numerical 
routine associated to QuaLiKiz has to be parallelized. 

The very first interest of running QuaLiKiz within an integrated transport solver is to 
analyze a tokamak discharge at a given time corresponding to quasi steady-state conditions. 
The main goal is the following: the radial profile of the effective heat and particle diffusiv- 
ities, as well of the temperatures and the density obtained from the interpretative analysis, 
should be matched by the predictive transport simulation using QuaLiKiz, provided that the 
heat and particle sources are constant in time. It is expected that, by the effect of the Qua- 
LiKiz transport coefficients, the temperature and density profiles and their relative gradients 
undergo to a time evolution, which should converge to the experimental values within the 
experimental uncertainty. More in particular, due to the absence of particle fueling in the 
center of the plasma, the steady-state condition corresponds to a null particle fiow in the 
core. 

A standard Tore Supra discharge, TS39596, already discussed in Chapter [3] and whose 
experimental parameters are summarized in Tables |4.1| and |4.2[ has been chosen as a first 
example to be analyzed with QuaLiKiz. The simulation has been run keeping a fixed density 
profile, while evolving both the ion and electron temperatures retaining the real Z^ff — 1.6. 
The transport coefficients predicted by QuaLiKiz are applied in the normalized radial domain 
< /9 < 0.8; for the outer plasma region corresponding to p > 0.8, an arbitrary profile of 
diffusivity is assumed, providing suitable boundary conditions for the temperature profiles. 
The results presented in Figs. |5.6| and |5.7| show the comparison between the Tg and 7^ pro- 
files obtained by the QuaLiKiz simulation and the measurements from the diagnostics. The 
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CRONOS simulation has been run for a time t Ki te-, where physically te ~ 200 ms, such 
that these temperature profiles reproduce a steady state on a macroscopic transport time 
scale. A reasonable agreement between the QuaLiKiz predictions and the measurements is 



obtained for both Tg and T^. A careful examination of the ion temperature profile of Fig. 5.7 



reveals that the boundary condition for Xi assumed at /? = 0.8 is not optimal, as highlighted 
by the comparison of the reconstructed profile with the Charge Exchange data. This can 
warn about the actual sensitivity of the prediction with respect to the boundary conditions. 
On the other hand, the peaking of the temperature profiles in the inner region {p < 0.3) is 
clarified when looking at the heat diffusivities calculated by the quasi-linear model and shown 



in Fig. 5.8 For p < 0.3 in fact, the model finds no turbulent transport since the local plasma 
parameters are below the linear mode thresholds: only the neoclassical contributions are left, 
hence a barrier-like effect appears in the temperature profiles. This point can be regarded as 
a typical consequence of the hypothesis of local transport. This feature is commonly shared 
among all the local transport simulations and, more importantly, it is not linked to the quasi- 
linear approximation, since a nonlinear local approach would lead to a similar result (see for 
example Ref. A solution to this general problem is beyond the scope this thesis work 

and represents one of the important challenges that have to be addressed to progress in the 
understanding and the prediction of the tokamak turbulent transport. 




Figure 5.6: First example of the application of QuaLiKiz integrated into the CRONOS code 
for the predictive analysis of the discharge TS39596 a,t t = 7.668 s. In this case only Tg 
and Ti are evolved using the experimental value of Z^ff = 1.6. The simulation achieves a 
steady state solution on a time At w Ite- This plot shows the radial profiles predicted by 
QuaLiKiz compared to the experimental data from both Electron Cyclotron Emission and 
Thomson Scattering diagnostics. 
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Figure 5.7: Radial Ti profiles predicted by QuaLiKiz compared to the experimental data 
from the Charge Exchange diagnostic. 




P 

Figure 5.8: Radial profiles of the ion and electron energy diffusivities predicted by QuaLiKiz. 



Chapter 6 

Conclusion 



The development of a gyro-kinetic quasi-linear transport model for tokamak plasmas, ul- 
timately designed to provide physically comprehensive predictions of the thermodynamic 
relevant quantities, is a task that required tight links among theoretical, experimental and 
numerical studies. The framework of QuaLiKiz, which operates a reduction of complexity 
on the nonlinear self-organizing plasma dynamics, allows in fact multiple validations of the 
current understanding of the tokamak micro-turbulence. 

The main outcomes of this thesis derive from the fundamental steps involved by the for- 
mulation of the quasi-linear transport model, namely; (1) the verification of the quasi-linear 
response against the numerically computed nonlinear one, (2) the improvement of the sat- 
uration model through an accurate validation of the nonlinear codes against the turbulence 
measurements, (3) the integration of the quasi-linear model within an integrated transport 
solver. 

The kinetic description of the plasma response is believed to be crucial to correctly cap- 
ture the essential physics at play in the tokamak turbulent transport. The TijT^ dependence 
of the ion ITG and electron TEM linear instability threshold is chosen as relevant example. 
Both the analytical calculations and the gyro-kinetic simulations recognize, for the first time 
in the case of the electron modes, that the temperatures ratio dependence of these thresholds 
is intrinsically linked with the kinetic wave-particle resonance. Hence, the gyro-kinetic for- 
mulation is preferred in QuaLiKiz, with respect to an arbitrary closure of a fluid or gyro-fluid 
model. 

The system of equations that defines the quasi-linear expectations for the turbulent fiuxes 

of energy and particle is presented, as it is implemented in QuaLiKiz. At the level of this 
formulation, two possible resonance broadenings are identified. The first one is related to the 
kinetic wave-particle resonance and its finite broadening is necessary to fulfill the ambipolar- 
ity of the particle fluxes. The second one is linked to an intrinsic frequency spectral shape 
of the fluctuation potential around the real frequency of the unstable mode. The validity of 
the hypothesis of the quasi-linear response is for the first time systematically investigated 
against fully nonlinear gyrokinetic simulations, apart from any choice on the structure of 
the saturated potential. Three different levels are examined. Firstly, the comparison of 
the characteristic turbulence times allows to derive an estimation of a Kubo number. Sec- 
ondly, the phase relations between the transported quantities and the fluctuating potential 
are studied in the linear as well in the nonlinear phase. Finally, the overall quasi-linear over 
nonlinear ratio for the transport of both the energy and particle channels is analyzed. The 
main outcome is that, with an appropriate re-normalization, the quasi-linear approximation 
is able of reasonably recovering the nonlinear fluxes over a wide range of plasma parameters. 
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within a mean rms error of 13%. This result surprisingly extends also to conditions where the 
high instability drive would have suggested strong turbulence conditions that could not be 
captured by the quasi-linear physics. Nevertheless, the quasi-linear particle flux can deviate 
from the nonlinear expccations when dealing with strong inward flows driven by ion modes. 
Fortunately, such conditions very rarely appear in the experiments. 

The model for the fluctuating saturated potential assumed in QuaLiKiz results from the 
validation of the numerical and theoretical predictions against the measurements. The non- 
linear gyrokinetic simulations, using the GYRO code, are here shown for the first time to 
quantitatively reproduce complementary measurements from a standard Tore Supra ohmic 
discharge, namely: the radial profiles of the heat diffusivity Xcft and of the rms 5n/n, as 
well the local k0 and kr spectra. When properly reconstructing the instrumental response 
of the diagnostic system, the simulations help moreover to revise the apparent experimental 
turbulence anisotropy at scales comparable to the ion gyroradius. Also, the similarity of the 
fiuctuation kg power spectra between the measurements and both local and global, using the 
GYSELA code, gyrokinetic simulations, suggests a rather general character of the tokamak 
turbulence wave-number spectrum. At the scales of the order of the; ion gyroradius, the 
power law scaling |5n/c|^ ~ |(5(/>fc|^ oc k~^, found on many other devices, is also recovered by 
an analytical shell model for the drift- wave turbulence. Simple theoretical considerations are 
again successfiil when applied to the experimental and the numerical frequency spectra of 
the fluctuations. An intermediate statistics between diffusive and convective behaviors for 
the particle displacements in the turbulent field, is found to agree with the experiments and 
the nonlinear simulations on both the frequency spectral shape and the kg dependence of the 
spectral broadenings. 

Finally, using only one renormalization constant, the total quasi-linear fluxes predicted 
by QuaLiKiz arc compared the a wide mimbcr of nonlinear gyrokinetic simulations. When 
coupling the choices for the saturated potential with the quasi-linear response, the QuaLiKiz 
fluxes are shown to agree with the nonlinear predictions for the energy transport in the ion and 
electron channels, as well as for particle fluxes for a wide range of tokamak relevant plasma 
parameters. QuaLiKiz is now coupled to the integrated transport platform CRONOS; its 
flrst applications to the experiment are encouraging but demand further investigation from 
both the quasi-linear and nonlinear modeling. 

A number of challenging issues still remains to be tackled. In fact, as it has been demon- 
strated in this thesis work, developing a reduced transport model allows to explore a very 
wide range of theoretical, experimental and numerical issues of interest for the nuclear fusion 
research. A first priority is certainly linked to further assess the effectiveness of QuaLiKiz in 
predicting the tokamak confinement properties, solving the time dependent transport prob- 
lem within CRONOS. A partial list of other potential topics of study and improvement are: i) 
The choices for the saturated potential deserve additional comparisons with nonlinear simu- 
lations and experimental measurements dealing with different plasma scenarios, ii) The way 
of accounting for the subdominant unstable modes in the quasi-linear formulation can be 
refined, iii) The domain in which the quasi-linear approximation fails (particle flux, marginal 
conditions, strong ITG turbulence, onset of zonal/large scale flows, etc.) should be more 
precisely understood, especially in view of the experimental plasma conditions, iv) Further 
improvements to the quasi-linear model could be addressed, namely accounting for the ExB 
shear stabilization effects and for shaped plasma geometries, v) The quasi-linear model can 
be extended to deal with the transport of momentum, in order to perform fully predictive 
studies, vi) The numerical optimization within the CRONOS transport solver will contribute 
to provide a fast and reliable tool for tokamak plasma studies. 
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The Tj/Te dependence of linear 
ITG-TEM thresholds 



The fluid approach 



Advanced fluid theories have been widely used for deriving analytical threshold expres- 
sions of electrostatic unstable modes. A significant reference is the Weiland reactive two- 
fluids model |105L I104j , which accounts for the influence and the interactions of convection, 
compression, and thermalization of the plasma species. Within the Weiland fluid approach, 
equations for ion and electron continuity, trapped electron and ion energy are considered. 
The hierarchy of fluid equations has necessarily to be truncated by a closure, which actually 
represents the wave-particle resonance. The Weiland model adopts the so called RighiLeduc 
or diamagnetic heat flow closure following by that of Braginskii |14) . In the first approxima- 
tion, the parallel ion motion can be neglected; this is reasonable if the fastest growing mode 
fulfills (fcpi)^ « 0.1 [T04|, where k is the mode wave vector and pi is the ion Larmor radius. 
The electron density perturbation drif. is written as: 



net 



(1 



It) 



(A.l) 



where index t and p stand, respectively, for trapped and passing particles. Passing electrons 
are allowed to reach a Boltzmann distribution, i.e., they are assumed adiabatic according 
to (Jnep/riep — e.54>/Tf,. These assumptions are coupled with the quasi-neutrality condition 
(assuming a single hydrogenous ion species): 

5ni = 5nep + Sn^t (A. 2) 

The linear part of the fluid equations leads to the following dispersion relation for a two- 
stream instability |1051 [79l I104j , where finite Larmor radius FLR effects have been neglected: 
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with the resonant denominator for the species j in the form: 
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(A.3) 



(A.4) 
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The Ti/Te dependence of linear ITG-TEM thresholds 



A first unstable branch is achieved with the condition Ni < that corresponds to a mode 
propagating in the ion magnetic drift direction (ITG modes), while for Ni > N,, the mode is 
rotating in the electron drift one (TE modes). 

Within the Weiland model these modes are thought to decouple when Ni <^ N^ (pure 
ITG mode) or Ni ^ N^ (pure TEM); for each of these cases, the full dispersion relation splits 
into two quadratic equations. Imposing a null imaginary part to the solution u) = ujj, + ij 
provides the instability threshold condition. In case of ITG modes this procedure leads to 
the following threshold expressed: 
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(A.5) 



where the main Ti/T^ proportionality is carried by the last term. In the advanced fluid 
Weiland model, this ITG Ti/T^ threshold dependence is carried by the choice on the resonant 



denominator, restoring the same feature highlighted by the kinetic approach in 2.2.4[ A plot 
for this threshold behavior versus R/Ln is presented in Fig. A.l for different ratios Ti/Te- For 



most common values of normalized density gradients, i.e. R/Ln < 5, the ITG threshold raises 
when increasing Ti/T^. Conversely, for R/Ln > 5, the opposite T^/Ti scaling is predicted 
by Eq. (A.5). Following an analogous procedure, the Weiland fluid model provides an 
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Figure A.l: ITG threshold vs R/Ln derived according to the Weiland model (Eq. (A.5)) for 
different values of the ratio Ti /T^ . 



analytical expression for the R/Ltc. TEM linear threshold: 
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(A.6) 



No term dependent on the temperature ratio is present within this formulation, contrary to 



the ITG case Eq. (A.5 1 



The Vrn/n TEM threshold in the limit of zero temperature gradients 



The very simplified case of V^T^ = V rTi = is addressed here, since trapped electron 
modes can be destabilized by the only presence of density gradients. Even if these are not 
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relevant conditions for tokamak plasmas, the present analysis is useful in order to enlarge 
predictive capabilities on wider parametric space. 

Dealing with zero temperature gradients implies that the linear kinetic dispersion relation 
(2.79) can be greatly simplified, since — uj^^ = 0. In this case it is not necessary to 



isolate the imaginary resonant contribution, which would lead to w* = (see Eq. (2.81)), 
where no energy dependence is present. The crude second order fluid expansion based on the 
condition w^ie/w «C 1 is then allowed here for a mode in the electron diamagnetic direction; 
the non-negligible ion response has been treated through a second order fluid expansion too. 
The critical R/Ln value above which unstable modes set in can only be numerically derived 
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Figure A. 2: T^/Ti dependence of the R/Ln TEM instability threshold calculated with 
Kinezero at R/Lti — R/Ltc — and at fixed Ti and T^. Analytical predictions follow- 
ing from a fluid expansion are also plotted. 



from the quadratic dispersion relation: 



(A.7) 



Ti \ ^ ^ ^ ^ / £ \ 

This Tf,/Ti dependence of R/ Ln,th has been compared with Kinezero simulations in Fig. 



A. 2 choosing the set of plasma parameters defined by the Table 2.1 and R/Lti = R/Ltc = 0. 
The fluid expansion foresees a weak decrease of the R/Ln TEM threshold for higher T^/Ti. 
On the other hand, linear gyrokinetic simulations identify a temperature ratio impact in 
agreement with analytical predictions only for low values of T^/Ti. When considering condi- 
tions of VrTe — '^rTi = 0, the crude fluid expansion based on nonresonant effects can then 
be regarded as an adequate treatment, capable of qualitatively reproducing the numerical 
features found by linear gyrokinetic simulations. Moreover, in the absence of temperature 
gradients, the ion response contribution has a weak impact of the ratio T^/Ti on the R/Ln 
TEM threshold. 
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The Ti/Te dependence of linear ITG-TEM thresholds 



Derivation of a consistent electron (ion) response within a kinetic approach 



The linear kinetic dispersion relation written in Eq. (2.76) can describe both ion and 
trapped electron modes. The analytical approach here presented refers to a single unstable 
branch in the ion or in the electron diamagnetic direction. The key point is evaluating the 
response of the opposite species with respect to the one giving the modes direction. In this 
appendix we examine in detail the case of electron modes; ion modes relations will not be 
here explicitly rewritten since they can be analogously obtained simply reversing the role of 
the species and accounting for the passing ions. 



As already explained in paragraph 2.2.4[ imposing a null imaginary part coming from the 
kinetic resonance provides fundamental conditions for the modes frequency (Eq. (2.801 for 



electron modes, Eq. (2.831 for ion ones) that can be rewritten for electron modes as: 
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The term proportional to R/Ln can be eventually set to zero when considering strict flat 
density limit a;*^ = 0. It is worth noting that the role of the density gradient inside expression 



(A. 8 1 limits the validity of this kind of approach entirely based on the kinetic resonance, 
because the modes frequency is not allowed to change sign. Kinetic effects will then be 
dominant on the threshold behavior until the conditions 
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th 



are fulfilled; beyond these limits the hypothesis of modes frequency close to the resonance is 
not valid and this procedure is not applicable anymore. 



Since the modes frequency is assumed to obey to Eq. ( A.8 1, it is possible to evaluate the 



ion response at the same frequency; ion resonant effects are in fact excluded because of the 



opposite sign of uj. On the other hand, the real part of Eq. (2.761 leads to 
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(A.10) 



The parameter H represents in fact the ion response in both the adiabatic and non-adiabatic 
components; it can be written as: 



H = 1 



+ 



= 1 - 



(A.11) 



The term proportional to R/Lti can be eventually set to zero in the absence of ion temper- 
ature gradients, w^- = 0. Physically, Eq. (A.11) expresses the ion response in its adiabatic 



and non-adiabatic contributions driven by ion temperature and density gradients. The only 
unknown is the normalized frequency uj/ujoTe', within our approach, the latter one is assumed 
to obey to Eq. (|A.9|). 

, Eq. 

which has to be evaluated according to Eq. (A.11); the un- 



Apart from R/Ln, Eq. (A.9) appears depending on the electron mode threshold itself 
through the term uj^^/ujDTe 



known threshold value and the ion response result then intrinsically linked one to each other. 
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Both Eqs. (A. 9) and ( A.10[ ) can be substituted in Eq. (A. 11 1, obtaining an implicit integral 
expression for H. Auto-consistent solutions for H can be finally numerically found, depending 
on the parameters T^/Te, i?/L„, R/Lti, and ft externally imposed. 

Within the already mentioned limits, this kind of procedure is then rigorously valid; the 
only approximations regard having neglected FLR effects and the simplification of the A- 
integration. 

Self-consistent calculations of Eqs. (A.9|,( A.10l,( A.ll I, considering w^,- = a;*^ = 0, have 
been carried out for deriving the Tf./Ti dependence of the pure TEM threshold in strict flat 
density limit: the result is shown in Fi g. |2.6[ Retaining a;pj,a;*g 7^ has instead lead to the 
analytical expectations shown in Fig. 2.7 An analogous set of coupled equations as Eqs. 
(A.9|,( A.10l,( A.ll I can be easily written when considering ion modes. With a;^g,a;*j ^ 
this has lead to the results shown in Fig. |2.9| 
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Brief notes on the nonlinear 
gyrokinetic codes 



About the GYRO code 

GYRO is a nonlinear tokamak micro-turbulence code. Developed at General Atomics 
since 1999, GYRO is presently one of the most advanced and comprehensive tool available 
for investigating the tokamak turbulence transport. GYRO uses a fixed Eulerian grid to 
solve the 5-D gyrokinetic-Maxwell equations using a 6f scheme. Despite the complexity of 
the nonlinear problem which requires massive parallel computing, the main feature of this 
code relies on its flexible operation, with the capability to treat: 



A local (otherwise called flux-tube) or global radial domain, in a full or partial torus. 
Generally shaped or simple circular s — a magnetic geometry. 

• Full nonlinear gyrokinetic treatment of ions, electrons and impurity species (both 
trapped and passing domains). 

• Electrostatic or electromagnetic fluctuations. 

• Electron-ion and ion-ion pitch-angle scattering operators. 

The right-handed, field-aligned coordinate system {ip,9,Q together with the Clebsch repre- 
sentation [5S] for the magnetic field is used. This has been briefly introduced in paragraph 



2.1.1 see in particular Eq. (2.4|, in the simplified case of concentric circular flux surfaces. 
Keeping the magnetic field representation B — V-q x V'0, more generally, the angle C can be 
written in terms of the toroidal angle as: 

C^^ + vii^.O) (B.l) 

In these coordinates, the formal Jacobian is : 

VV' X V0 • VC VV" X • \dr ) ^ ' ' 



The advantage of this representation is the capability of treating a general shaped magnetic 
equilibrium, solution of the Grad-Shafranov equation, and not only circular flux surfaces. 
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Brief notes on the nonlinear gyrokinetic codes 



Recalling the usual expression for the magnetic field from Eq. (2.1 ), v [ip, 0) is expressed by 
an integral relation: 

i^i^,e) = -Ii^) J^\V^fd0 (B.3) 
Jo 

While for a general Grad-Shafranov equilibrium, v {')p,9) can be only numerically solved, in 
the simplified case of concentric circular s — a geometry, we recover the more usual result 
I' {ifj, 9) = —q (ip) 9 and therefore C, — ip — q (ip) 9. 

Finally in this framework, the following flux-surface operator acting on the general quantity 
z is defined: 

dC I d9J^z 

= (B.4) 



/ dC / d9J^ 
Jo Jo 



The Sf expansion refers to the single particle distribution function /s, which is written 
as a sum of an equilibrium part and fluctuating terms: 

fs (x, E, t) = fos (x, E) + Sf, (x, E, M, t) (B.5) 

where x = R + p is the particle position, p is the gyroradius vector and R is the guiding- 
center position. The equilibrium function /os is assumed to be a Maxwellian. 
In GYRO, all the perturbed quantities are expended as Fourier series in C. For example, the 
normalized electrostatic potential Sep = e5(j)/Tg is written as: 

50(r,0,C)= £ <50„(r,0)e™^ n = j An (B.6) 

Contrarily to other nonlinear codes, which use different discretization schemes, GYRO op- 
erates through a direct discretization of the quantities z {r,9). From the numerical point of 
view, GYRO always solves the discrete linear or nonlinear gyrokinetic and Poisson/ Ampere 
equations in {n,r,9), plus the time advance. One of the perpendicular coordinates is then 
the radial one r (otherwise referred as x) ; the other perpendicular direction (usually referred 
as y) is not a coordinate in GYRO, which uses instead the Fourier index n, conjugate to the 
Clebsch angle C ~ V + (''■ (^)- 

Considering a physical function represented by the real field z{r,9,(^), a 9 periodicity 
condition can be written as: 

z{r,0,ip + iy (r, 0)) = z (r, 2TT,ip + i' (r, 2n)) (B.7) 

So, even if the physical field z is 27r-periodic in 9, the Fourier coefficients of the expansion 



(B.6) are not periodic, while they satisfy the phase condition: 

z„(r,0) = e2™«Wz„(r,2^) (B.8) 
Moreover, since z is real, the Fourier coefficients satisfy the relation — The spectral 



form of Eq. (B.6) is then 27r/An-periodic in ( (or in at fixed (r, 9). A direct consequence 



of the representation (B.6) allows to easily map the physical quantities at the outboard 



midplane in terms of the coordinates (r, ip) : 

N„-l 

E 



z{r,9 = 0,ip,t) = y z„ (r, 61 = 0, e-*""^ (B.9) 
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The perturbed distribution function Sfg solved by GYRO is given by (Gaussian CGS 
units are here used): 



Sfs (x,£;,/i,i) 



fos 



. (x, t) - g'5(j) (R, t) + ^e^'Mii (R, t) + h' (R, E, /i, t) (B.IO) 



Using the following notation 
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(B.ll) 

(B.12) 

(B.13) 
(B.14) 

(B.15) 



where G is the gyroaverage operator, the gyrokinetic equation used in GYRO for a generic 
species s can be synthetically written as: 



dhs 
dt 



wus ■yhs^ + inQs^*sQ''Us^ = -v\\sV\\Hs - vds ■ ^Hs + C [ h, 

1 2 



(B.16) 



The driver perturbed field Q'^Us is given by the Poisson- Ampere equation: the perturbed 
charge and current densities determining Q^Us enter as velocity space integrals over hs, 
which is proportional to nps — ^Os/'^oej the density of the species relative to the electron 
density. In the gyrokinetic Eq. (B.16), the terms 1 and 2 are respectively the nonlinear and 
linear generalized E x B drift motions, including the magnetic fluctuations. The terms 3 and 
4 represent instead the parallel and curvature drift motions, while the term 5 accounts for 
the collisions. 

Finally, the particle and energy flux for each species are deflned as: 



Ts ir)=^J dv^^fs (x) (^^h X VC/^ • Vr 
Qs i"")^^/ dv^E6fs (x) (J^h X VC/^ • Vr 



(B.17) 
(B.18) 



Tracer and quasi-linear transport 



A species is here qualifled as tracer when fios ^ 1- The following statements are derived 
and can be verified by GYRO simulations: 

1. The tracer species have a negligible contribution to the fields through the Poisson- 
Ampere equation; hence, the background plasma turbulence and transport is unaffected 
by the presence of the tracers. 

2. If the tracer species have identical mass and charge to the main plasma species, and 
all the terms of the gyrokinetic equation (B.16) for tracers are kept, then the tracer 
particle and energy fiuxes are identical to those of the main plasma species. 
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3. The tracer particle and energy fluxes, unlike the main plasma species ones, are indepen- 
dent of the tracer gradients \/Ln,tr and l/LT,tr as long as the tracer rjtr = Ln,tr/LT,tr 
is held fixed. The tracer transport follows in fact a Pick's law, i.e. linear with the 
gradients only at fixed rjtr- 



The full frequency spectrum quasi-linear approach (fQL) described in the paragraph 3.2.3 
refers to GYRO tracer (test-particle) simulations. Apart from the main plasma species, ion 
and electrons, two additional tracers are retained, with identical masses and charges to the 
main ions and electrons. Conversely, these tracers are characterized by negligible densities, i.e. 
noe,tr = noi,tr = 10~^. Morcovcr, in order to deal with quasi-linear test-particle transport. 



the nonlinear term 1 in the gyrokinetic equation (B.16 1 governing both the tracers is dropped 



As highlighted in the paragraph |3.2.3[ this method should allow to obtain the best quasi-linear 
transport estimate. The main advantages of the present fQL approach are in fact: (1) the 
whole structure of the linear modes (and not only the linear most unstable and ballooning 
mode) is retained within the calculation, (2) the nonlinear saturation of the fluctuating 
potential is completely and self-consistently imposed by the main plasma species. On the 
other hand, the main drawback of the method is that the particle flux ambipolarity is not 
respected anymore. 

Apart from the numerical tracers experiments, an alternative approach can be used to 




400 600 
time [a/c ] 



1000 



Figure B.l: GYRO simulation on the GA-ITG-TEM standard case using the quasi- linear 
test-particle method: the time evolution of the energy diffusivities is shown for both the 
main plasma and the quasi-linear tracer species. 



investigate the quasi-linear transport. For a standard ions and electrons plasmas, recalling 
the gyrokinetic equation Eq. (B.16), only the linear E x B drive expressed by the term 1 
can be artificially suppressed in a nonlinear simulation for a single toroidal wave-number. 
Consequently, there will be no linearly unstable mode in the deleted wave number but still 
the nonlinear dynamics. Fig. |B.2| reports the results of this exercise with the GA-ITG-TEM 
standard case, where the linear drive has been suppressed for the wave-numbers corresponding 
to kePs — 0.1 and kgps — 0.6. As shown in the figure, some transport appears also where 
the linear drive is deleted, due to purely nonlinear couplings in the k space. It is worth 
noting that in this case, the more relevant purely nonlinear contribution in Fig. B.2 refers to 
the particle transport, i.e. the channel where the quasi-linear response appears more feeble 
(see paragraph 3.2.3). Analogous demonstration of an high-fc transport nonlinearly driven 
by low-fc scales has been investigated in Ref. [100] . 



Ill 




0.4 



0.4 



Figure B.2: GA-ITG-TEM standard case ion energy (a), electron energy (b), and particle 
(c) contribution per mode diffusivity versus wave number with the linear E x B drive at 
keps = 0.10 and kgps — 0.60 deleted hence linearly stable. Purely nonlinearly induced 
transport in highlighted in red. The total transport is somewhat reduced. 



Coordinates systems 

As previously mentioned, GYRO use a field aligned coordinate systems and it does not em- 
ploy an expansion in poloidal harmonics. Since several nonlinear tokamak micro-turbulence 
codes adopt different coordinates, it is of interest to detail some basic relations useful in 
order to make comparisons between different representations. In particular, here we refer to 
the coordinate system used in the GYSELA code. The latter one operates in fact with the 
coordinates (r, 0, ip) (respectively the radial coordinate, the poloidal and the toroidal angles): 
conversely to GYRO, the fields are expanded in Fourier series in both 9 and Lp. 

In GYSELA, the expansion for a generic field z is written as: 

.«^(r,^,V.) = ^z„«^e»("-+™«) (B.19) 

n,m 

which has to be compared to the analogous expression for GYRO. The superscripts GR and 
GS are here introduced to label GYRO and GYSELA notations respectively. The following 
relations link the GYSELA representation with the GYRO one within a simplified circular 
magnetic s — a equilibrium, where the relation C = + (r, 0) reduces to the analytically 
tractable form ( = (p — q (r) 9. A fundamental relation is then: 

J2 ^ (r, 9) e"< = (r) e<-^+"^<'^ {B.20) 

n n,m 

hence it follows: 

ir) = r ^z^^ (r, 9) e-I^+^^C^)!^ (B.22) 

which can be practically used to map one representation into another one. 

When dealing with quantities that are defined through the flux-surface average defined 
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by Eq. (B.4|, it appears that: 

Q^"" (r) - l[ f {r, e) (r, 9) (B.23) 
-f§^T. (r) (r) e<--y = ^ z^^ (r) vfj (r) (B.24) 

^ m.m' m 

The latter relation is useful when defining for example a flux contribution coming from a 
certain toroidal wave-number n. It also immediately leads to the definition of a n-dependent 
flux-surface averaged rms quantity, i.e. 



( ^ 

GR , 



Appendix C 



The case of quasi-linear pure 
ITG turbulence 



Even if realistic application of the transport models to tokamak plasmas demands that the 
non-adiabatic electron physics is retained, it is still relevant to test the validity of the quasi- 
linear approximation in the case of pure ion ITG turbulence, i.e. considering adiabatic 
electrons. This has been done again using the GYRO code in the local limit, adopting the 
same parameters and numerical resolution summarized in Table [XT] except for the electron 
response which is forced to be completely adiabatic. 

As previously done, the validity of the quasi-linear response is tested with GYRO through 
both the approaches mQL, i.e. retaining only the linear leading mode, and fQL, i.e. the 
full frequency spectrum approach (see paragraph 3.2.3). A wide scan of the normalized 
ion temperature gradient 1.5 < a/L^ < 9.0 is performed and the results are reported in 
Fig. |C.1| A clear breakdown of the mQL approximation is found for the ion energy flow 
driven in the case of pure ITG, going far from threshold to very high turbulence levels. Fig. 
C.l| shows that quasi-linear over nonlinear ratio (overage) computed according to the mQL 



approach increases up by a factor of 2.1 from the overage of 1.64 at the reference value 
of a/L^i = 3.0. Surprisingly, the quasi-linear over nonlinear ratio computed with the fQL 
approach does not exhibit a similar feature, while it stays reasonably constant across the 
whole scan. 

Relevant insight can be gained when studying the cross-phase relations, reported in Fig. 
|C.2[ comparing the phase angles of the nonlinear saturation regime and those of the linear 
most unstable mode. Still surprisingly, the linear cross-phases relative to the leading linear 
mode accurately track the nonlinear ones, even for the highest turbulence levels. The results 
of Fig. |C.2| are then a clear example highlighting that the information of the quasi-linear 
response is not completely carried by the cross-phases. In the case of pure ITG turbulence in 



fact, the failure of the mQL approximation reported in Fig. C.l is linked to the quasi-linear 
relative amplitudes \5Ei \ / \5(j)\, which originate an over-estimate of the total turbulent energy 
flux. A relevant question at this point is then: why the fQL approach does not lead to a 
similar failure? A possible explanation is here below proposed. 

The effects coming from the toroidicity are very strong for pure ITG turbulence. The guess 
is that a relevant physical mechanism linked to this point is missing in the mQL approach, 
which is instead properly retained with the fQL approximation. It has been already discussed 
that mQL deals with only the linear leading mode: more in particular this always refers to a 
ballooning normal mode, i.e. a mode centered around a fixed angle (this angle is usually. 
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Figure C.l: Quasi-linear over nonlinear overages across a ion temperature gradient scan based 
for pure ion ITG turbulence, using both the niQL and the fQL approaches. Xi are also given 
in terms of XGB units. 




-2 2 -2 2 -2 2 

angle [rad] angle [rad] angle [rad] 



Figure C.2: PDF of the nonlinear cross-phases and the linear cross-phase of the most unstable 
mode (white line): a) 6Ei — Scj) at a/Lri — 2.0, b) SEi — Scf) at a/LTi — 3.0 and c) SEi — Scj) at 
a/LTi — 9.0, from local GYRO simulations of pure ITG turbulence with adiabatic electrons. 
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but not necessarily, = 0). Normally, other less unstable modes at different ballooning 
angles are active and contributing to the nonlinear saturation. Previously, it has been shown 
how in the mQL approach, the neglect of these less ballooning and less unstable modes can 
account for the different quasi-linear over nonlinear ratios seen between the mQL and the 
fQL approximations for coupled ITG-TEM turbulence, as reported in Fig. |3.7| 

This argument can be made more explicit recalling the general expression for the quasi- 
linear turbulent flux ( 3.37[ ). A more general relation can in fact be written as: 



QL-flux cx QL-weight;. J ® Spectral-intensity ^.^ g^j (C.l) 
where also the convolution over the ballooning angles is retained. The mQL approach 



results from the Eq. (C.l) in the limit of both j — j and Oq — !■ Oq, meaning respectively 
the most unstable j mode centered at the ballooning angle. On the other hand, the fQL 
approximation correctly retains the quasi- linear convolution over j and ^o- 
Hence, a possible origin of the failure of the mQL approximation for pure ITG turbulence 
is that the quasi-linear weight is computed only at a fixed ballooning angle 9q. Since the 
nonlinear spectral intensity can not be separated into distinguishable contributions in j and 
Oq, this implies that the saturation potential spectrum is entirely applied to the single bal- 



looning mode do instead of the proper convolution of Eq. (C.l ). 

Finally it can be concluded that the validity of the quasi-linear approximation has been 
successfully tested in the case of pure ITG turbulence across a wide scan of the ion temper- 



ature gradient, using the fQL approach (Fig. C.l). Nevertheless it has to be stressed that 



this requires accounting for the whole spectrum of the 6q ballooning linear modes. 
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